## This library contains functions I normally use for Bayesian simulation ##########################################################################

######################## Index of functions####################################
# 1. AcceptRate        # Function calculate acceptance probability in MH for AR(2)
# 2. AcceptRateARp     # Function calculate acceptance probability in MH for AR(p)
# 3. AR2Cov            # Computing the covariance matrix of AR(2) process
# 4. BBDraw            # Sampling bi for the conditional likelihood algorithm
# 5. BetaTerm          # Computing terms for conditional density of beta for the conditional likelihood algorithm
# 6. bdiag             # Creating a diagonal block matrix
# 7. BiNorm1           # Computing the the mean and variance of conditional distibution of first element in a bivariate normal distribution
# 8. BiNorm2           # Computing the the mean and variance of conditional distibution of second element in a bivariate normal distribution     
# 9. BTnorm            # sampling from a truncated bivariate normal distribution
# 10. Cov.AR1          # computing the covariance matrix of AR(1) matrix
# 11. Cov.AR2          # computing the covariance matrix of AR(2) matrix    
# 12. Cov.ARp          # computing the covariance matrix of AR(p) matrix 
# 13. crosssum         #  sum up the term of X%*% Y for all the observations
# 14. CutPointUpdate   # sampling the cutpoints for ordinal model
# 15. DataAug          # conducting data augmentation in the conditional likelihood algorithm
# 16. GroupPred        # computing the mean of variables for group-level predictors
# 17. Index            # generating index of groups

# 18. myrwish          # robust version of random sample from wishart distribution by using sechol instead of chol

# 19. myriwish         # robust version of random sample from inverse wishart distribution by using sechol instead of chol

# 20. my.positive.making # robust version of making a positive definite matrix by taking care of complex eigen values

# 21. OrdinalAugment   # data augmentation for ordinal data
# 22. outerself        # outer product of a matrix
# 23. PhiTerm          # computing the proposal density of phi, outdated version
# 24. PhiTerm2         # computing the proposal density of phi, good version
# 25. PhiTermARp       # computing the proposal density of phi, AR(p)
# 26. PhiProp          # sampling proposal from the stationarity space: AR(2)
# 27. PhiProp          # sampling proposal from the stationarity space: AR(p)
# 28. Poly             # transforming data by using polynomial
# 29. random           # computing the interacting matrix for random effects
# 30. Random.Coef      # computing the random coefficients of beta (unit)
# 31. Random.Coef.gamma # computing the random coefficients of gamma (time)
# 32. Rescale          # rescaling data such that zero-mean and 1 variance
# 33. SimRandTime      # sampling ct in the conditional likelihood algorithm
# 34. ss               # computing the term in the exponential for the acceptance probability
# 35. sumtrans         # sum the term for all the observations --- for use of posterior variance (normal) update
# 36. trans2           # transform the data for phi
# 37. TRFirst          # data augmentation for the first observation: conditional likelihood version
# 38. TRSecond          # data augmentation for the second observation: conditional likelihood version
# 39. TRSecondLast      # data augmentation for the second last observation: conditional likelihood version
# 40. TRSecondl2        # data augmentation for the second last observation: conditional likelihood version
# 41. TRSecondl3        # data augmentation for the second last observation: conditional likelihood version
# 42 TRThird            # data augmentation for the third observation: conditional likelihood version

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

options(digits = 5) 
library(MCMCpack);  library(mvtnorm);  library(lattice); library(msm)
library(MASS); #library (XunLibrary1)
#library(agce);
library(Matrix); library(corpcor); library(accuracy)
library(QuantPsyc);  library(panel); library(kinship); library(bayesSurv)

AcceptRate <- function(rho1=phi1.draw, rho2= phi2.draw, orho1=phig[1], orho2=phig[2], Osig=Sig , unit=id1, AY=upystar, AX=X, AW=W, Abeta=betadrawn, Ab=ttn, AT=id2, Ac=ccin, Entry=entry, f=g){
  
    newphi<-rbind(rho1, rho2)
    oldphi<-rbind(orho1, orho2)
    pG<-matrix(c(rho1, 1, rho2, 0), ncol=2)
    ll<-II-kronecker(pG,pG)
    pd<-as.vector(pseudoinverse(ll)%*%as.vector(e%*%t(e)))
    newsig<-matrix(pd, nrow=2)
      nSig <- det(newsig)
      oSig <- det(Osig)
      fphi1 <- 0
      ofphi1 <-0
      fphi <- 0
      ofphi <-0
    
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
    T<-AT[roww]
    t<-min(T)
   if (length(roww)==1){
     fphi1 <- fphi1-1/2*log(newsig[1,1])-(AY[roww]-AX[roww,]%*%Abeta-AW[roww,]%*%Ab[n,]-Ac[t])^2/(2*newsig[1,1])
     ofphi1 <- ofphi1 -1/2*log(Osig[1,1])-(AY[roww]-AX[roww,]%*% Abeta-AW[roww,]%*%Ab[n,]-Ac[t])^2/(2*Osig[1,1])
     }else{
    P<-roww[c(1,2)]
    ystari<-AY[P]
    Xi<-AX[P,]
    Wi<-AW[P,]
    bbbn<-Ab[n,]
    ccn <- cbind(Ac[ c(t, t+1)])
    fphi <- fphi-1/2*log(nSig)+ss(ff=newsig, Yb=ystari, Xb=cbind(Xi, Wi, ccn), bb=c(Abeta, bbbn, 1))
    ofphi <- ofphi -1/2*log(oSig)+ss(ff=Osig, Yb=ystari, Xb=cbind(Xi, Wi, ccn), bb=c(Abeta, bbbn, 1))
       }
     }
   ppalpha<-fphi1+fphi-ofphi1-ofphi
   palpha<- exp(ppalpha)
    alpha <- min(palpha, 1)
      return(alpha)
  }

AcceptRate1 <- function(rho=phi.draw, orho=phig,precision=1, unit=id1, AY=upystar, AX=X, AW=W, Abeta=betadrawn, Ab=ttn,  Entry=entry){
  
   fphi <- 0
   ofphi <- 0
   
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
    T<-AT[roww]
    t<-min(T)
    P<-roww[1]
    ystari<-AY[P]
    Xi<- t(as.matrix(AX[P,]))
    Wi<- t(as.matrix(AW[P,]))
    bbbn <- as.matrix(Ab[n,])
    fphi <- fphi+1/2*log(precision*(1-rho^2))-precision*(1-rho^2)/2*(ystari-Xi%*%cbind(Abeta)-Wi%*%cbind(bbbn))^2    # check: why I used + ?
   ofphi <- ofphi+1/2*log(precision*(1-orho^2))-precision*(1-orho^2)/2*(ystari-Xi%*%cbind(Abeta)-Wi%*%cbind(bbbn))^2
  }
   ppalpha<-fphi-ofphi
   palpha<- exp(ppalpha)
   alpha <- min(palpha, 1)
      return(alpha)
  }

AcceptRate0 <- function(rho=phi.draw, orho=phig,precision=1, unit=id1, AY=upystar, AX=X, AW=W, AS=S, Abeta=betadrawn, Ab=ttn, AT=id2, Ac=ccin, Entry=entry){
  
   fphi <- 0
   ofphi <- 0
   
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
    T<-AT[roww]
    t<-min(T)
    P<-roww[1]
    ystari<-AY[P]
    Xi<- t(as.matrix(AX[P,]))
    Wi<- t(as.matrix(AW[P,]))
    Si <- t(as.matrix(AS[P,]))
    bbbn<- as.matrix(Ab[n,])
    ccn <- as.matrix(Ac[t,])
       rant <- Si%*%cbind(ccn)
    fphi <- fphi+1/2*log(precision*(1-rho^2))-precision*(1-rho^2)/2*(ystari-Xi%*%cbind(Abeta)-Wi%*%cbind(bbbn)-rant)^2 
   ofphi <- ofphi+1/2*log(precision*(1-orho^2))-precision*(1-orho^2)/2*(ystari-Xi%*%cbind(Abeta)-Wi%*%cbind(bbbn)-rant)^2
  }
   ppalpha<-fphi-ofphi
   palpha<- exp(ppalpha)
   alpha <- min(palpha, 1)
      return(alpha)
  }

AcceptRateARp <- function(nrho=phi.draw, orho=phig, Osig=Sig , unit=id1, AY=upystar, AX=X, AW=W, precision=1, Abeta=betadrawn, Ab=ttn, AT=id2, Ac=ccin, Entry=entry, f=g, lag=p){
  
    II <- diag(lag^2)
    hh <- rep(0, (lag-1))
    e <- c(1, hh)
    newphi<-nrho
    oldphi<-orho
    pG<-matrix(NA, ncol=lag, nrow=lag)
    pG[2:lag, 1:(lag-1)] <- diag(lag-1)
    pG[1,] <- newphi
    pG[-1, lag] <- 0
     
     ll<-II-kronecker(pG,pG)
     pd<-as.vector(solve(ll)%*%as.vector(e%*%t(e)))
   
    newsig<-1/precision* matrix(pd, nrow=lag)
    nSig <- det(newsig)
    oSig <- det(Osig)
    fphi1 <- 0
    ofphi1 <-0
    fphi <- 0
    ofphi <-0
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
    T<-AT[roww]
    t<-min(T)
    Ti <- length(roww)
    if (Ti==1){
     fphi <- fphi-1/2*log(newsig[1,1])-(AY[roww]-AX[roww,]%*%Abeta-AW[roww,]%*%Ab[n,]-Ac[t])^2/(2*newsig[1,1])
     ofphi <- ofphi -1/2*log(Osig[1,1])-(AY[roww]-AX[roww,]%*% Abeta-AW[roww,]%*%Ab[n,]-Ac[t])^2/(2*Osig[1,1])
      }else{
      Ji<-min(Ti, lag)
      P<-roww[c(1:Ji)]
      ystari<-AY[P]
      Xi<-AX[P,]
      Wi<-AW[P,]
      bbbn<-Ab[n,]
      ccn <- cbind(Ac[ c(t, t+lag-1)])
     nnSig <- newsig[1:Ji, 1:Ji]
     ooSig <- Osig[1:Ji, 1:Ji]
     dnSig <- det(nnSig)
     doSig <- det(ooSig)
     bbbn<-Ab[n,]
     ccn <- cbind(Ac[ c(t:(t+Ji-1))])
     fphi <- fphi-1/2*log(dnSig)+ss(ff=nnSig, Yb=ystari, Xb=cbind(Xi, Wi, ccn), bb=c(Abeta, bbbn, 1))
     ofphi <- ofphi -1/2*log(doSig)+ss(ff=ooSig, Yb=ystari, Xb=cbind(Xi, Wi, ccn), bb=c(Abeta, bbbn, 1))
 }
   }
   ppalpha<-fphi-ofphi
   palpha<- exp(ppalpha)
    alpha <- min(palpha, 1)
      return(alpha)
  }


AcceptRateARp1 <- function(nrho=phi.draw, orho=phig, Osig=Sig , unit=id1, AY=upystar, AX=X, AW=W, Abeta=betadrawn, Ab=ttn, Entry=entry, f=g, lag=p){

    II <- diag(lag^2)
    hh <- rep(0, (lag-1))
    e <- c(1, hh)
     newphi<-nrho
    oldphi<-orho
    pG<-matrix(NA, ncol=lag, nrow=lag)
    pG[2:lag, 1:(lag-1)] <- diag(lag-1)
    pG[1,] <- newphi
    pG[-1, lag] <- 0
     
    ll<-II-kronecker(pG,pG)
    pd<-as.vector(pseudoinverse(ll)%*%as.vector(e%*%t(e)))
    
    newsig<-matrix(pd, nrow=lag)
    nSig <- det(newsig)
    oSig <- det(Osig)
    fphi1 <- 0
    ofphi1 <-0
    fphi <- 0
    ofphi <-0
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
     Ti <- length(roww)
    if (Ti==1){
     fphi <- fphi+1/2*log(newsig[1,1])-(AY[roww]-AX[roww,]%*%Abeta-AW[roww,]%*%Ab[n,])^2/(2*newsig[1,1])
     ofphi <- ofphi +1/2*log(Osig[1,1])-(AY[roww]-AX[roww,]%*% Abeta-AW[roww,]%*%Ab[n,])^2/(2*Osig[1,1])
   }else{
     Ji<-min(Ti, lag)
      P<-roww[c(1:Ji)]
      ystari<-AY[P]
      Xi<-AX[P,]
      Wi<-AW[P,]
      bbbn<-Ab[n,]
      nnSig <- newsig[1:Ji, 1:Ji]
      ooSig <- Osig[1:Ji, 1:Ji]
      dnSig <- det(nnSig)
      doSig <- det(ooSig)
     bbbn<-Ab[n,]
     fphi <- fphi-1/2*log(dnSig)+ss(ff=nnSig, Yb=ystari, Xb=cbind(Xi, Wi), bb=c(Abeta, bbbn))
    ofphi <- ofphi -1/2*log(doSig)+ss(ff=ooSig, Yb=ystari, Xb=cbind(Xi, Wi), bb=c(Abeta, bbbn))
 }
   }
   ppalpha<-fphi-ofphi
   palpha<- exp(ppalpha)
    alpha <- min(palpha, 1)
      return(alpha)
  }

# Fix=Xi; Rand=Wi;  Ldep=ystar; FP=betan; TRP=ccin; unit=n; Lag=p; rrowss=rows; vvar=var; RRB=RB; DDi=Di; phig=phig; TSig=SSig; id2=id2;  Xstar=Xstar; Wstar=Wstar; cstar=cstar
AcceptRateC1 <- function(nrho=phi.draw, orho=phig, Osig=Sig , unit=id1, AY=upystar, AX=X, AW=W, precision=1, Abeta=betadrawn, Ab=ttn, Entry=entry, f=g, lag=p){
    II <- diag(lag^2)
    hh <- rep(0, (lag-1))
    e <- c(1, hh)
     newphi<-nrho
    oldphi<-orho
    pG<-matrix(NA, ncol=lag, nrow=lag)
    pG[2:lag, 1:(lag-1)] <- diag(lag-1)
    pG[1,] <- newphi
    pG[-1, lag] <- 0
     
    ll<-II-kronecker(pG,pG)
     pd<-as.vector(pseudoinverse(ll)%*%as.vector(e%*%t(e)))
   
    newsig<-1/precision* matrix(pd, nrow=lag)
    nSig <- det(newsig)
    oSig <- det(Osig)
    fphi1 <- 0
    ofphi1 <-0
    fphi <- 0
    ofphi <-0
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
     Ti <- length(roww)
    if (Ti==1){
     fphi <- fphi-1/2*log(newsig[1,1])-(AY[roww]-AX[roww,]%*%Abeta-AW[roww,]%*%Ab[n,])^2/(2*newsig[1,1])
     ofphi <- ofphi -1/2*log(Osig[1,1])-(AY[roww]-AX[roww,]%*% Abeta-AW[roww,]%*%Ab[n,])^2/(2*Osig[1,1])
   }else{
     Ji<-min(Ti, lag)
    P<-roww[c(1:Ji)]
   ystari<-AY[P]
   Xi<-as.matrix(AX[P,])
    Wi<-as.matrix(AW[P,])
   bbbn<-as.matrix(Ab[n,])
     nnSig <- newsig[1:Ji, 1:Ji]
     ooSig <- Osig[1:Ji, 1:Ji]
     dnSig <- det(nnSig)
    doSig <- det(ooSig)
     bbbn<-as.matrix(Ab[n,])
    fphi <- fphi-1/2*log(dnSig)+ss(ff=nnSig, Yb=ystari, Xb=cbind(Xi, Wi), bb=c(Abeta, bbbn))
    ofphi <- ofphi -1/2*log(doSig)+ss(ff=ooSig, Yb=ystari, Xb=cbind(Xi, Wi), bb=c(Abeta, bbbn))
 }
   }
   ppalpha<-fphi-ofphi
   palpha<- exp(ppalpha)
    alpha <- min(palpha, 1)
      return(alpha)
  }

AcceptRateARpT <- function(nrho=phi.draw, orho=phig, Osig=Sig , unit=id1, AY=upystar, AX=X, AW=W, AS=S, precision=1, Abeta=betadrawn, Ab=ttn, AT=id2, Ac=ccin, Entry=entry, f=g, lag=p){
  
    II <- diag(lag^2)
    hh <- rep(0, (lag-1))
    e <- c(1, hh)
    newphi<-nrho
    oldphi<-orho
    pG<-matrix(NA, ncol=lag, nrow=lag)
    pG[2:lag, 1:(lag-1)] <- diag(lag-1)
    pG[1,] <- newphi
    pG[-1, lag] <- 0
     
     ll<-II-kronecker(pG,pG)
     pd<-as.vector(solve(ll)%*%as.vector(e%*%t(e)))
   
    newsig<-1/precision* matrix(pd, nrow=lag)
    nSig <- det(newsig)
    oSig <- det(Osig)
    fphi1 <- 0
    ofphi1 <-0
    fphi <- 0
    ofphi <-0
   for (n in 1: max(unit)){
    roww <- Entry[which(unit==n)]
    T<-AT[roww]
    t<-min(T)
    Ti <- length(roww)
    if (Ti==1){
     fphi <- fphi+1/2*log(newsig[1,1])-(AY[roww]-AX[roww,]%*%Abeta-AW[roww,]%*%Ab[n,]-AS[rowww,]%*%timerandom[t,])^2/(2*newsig[1,1])
     ofphi <- ofphi + 1/2*log(Osig[1,1])-(AY[roww]-AX[roww,]%*% Abeta-AW[roww,]%*%Ab[n,]-AS[rowww,]%*%timerandom[t,])^2/(2*Osig[1,1])
      }else{
      Ji<-min(Ti, lag)
      P<-roww[c(1:Ji)]
      ystari<-AY[P]
      Xi<- as.matrix(AX[P,])
      Wi<- as.matrix(AW[P,])
      Si <- as.matrix(AS[P,])
      bbbn<-as.matrix(Ab[n,])
     nnSig <- newsig[1:Ji, 1:Ji]
     ooSig <- Osig[1:Ji, 1:Ji]
     dnSig <- det(nnSig)
     doSig <- det(ooSig)
     bbbn<-as.matrix(Ab[n,])
     ccn <- as.matrix(cbind(Ac[ c(t:(t+Ji-1)),]))
      rant <- c()
      for (k in 1:Ji){
        rant[k] <- Si[k,]%*%cbind(ccn[k,])
      }
     fphi <- fphi-1/2*log(dnSig)+ss(ff=nnSig, Yb=ystari, Xb=cbind(Xi, Wi, rant), bb=c(Abeta, bbbn, 1))
     ofphi <- ofphi -1/2*log(doSig)+ss(ff=ooSig, Yb=ystari, Xb=cbind(Xi, Wi, rant), bb=c(Abeta, bbbn, 1))
 }
   }
   ppalpha<-fphi-ofphi
   palpha<- exp(ppalpha)
    alpha <- min(palpha, 1)
    return(alpha)
  }

BBDraw<-function(Fix=Xi, Rand=Wi,  Ldep=ystar, FP=betan, TRP=ccin, unit=n, Lag=p, rrowss=rows, vvar=var, RRB=RB, DDi=Di, phig=phig, TSig=SSig, id2=id2, Xstar=Xstar, Wstar=Wstar, cstar=cstar){
     b.drawn<-c()
     T<-id2[rrowss]
      t<-min(T)
    if (length(rrowss)==1){
        b1 <- pseudoinverse(1/vvar*(cbind(Rand) %*% rbind(Rand))+pseudoinverse(DDi))
         b1 <- my.positive.making(b1)
         b1 <- matrix(as.real(b1), nrow=RRB)
        bmean1 <- b1%*%(1/vvar*cbind(Rand)%*%(Ldep-Fix%*% cbind(FP)-TRP[t]))
          b.drawn<-mvrnorm(1, bmean1, b1)
        }else{
          if (length(rrowss)==2){
            XP<-rbind(Fix[1,], Fix[2,])
            WP<-rbind(Rand[1,], Rand[2,])
            YP<-rbind(Ldep[1], Ldep[2])
     b1<-pseudoinverse((t(WP)%*%TSig%*%(WP)+pseudoinverse(DDi)))
     b1 <- make.positive.definite(b1)
     b1 <- matrix(as.real(b1), nrow=RRB)
     bmean1<-t(b1%*%(t(WP)%*%TSig%*%(YP-XP%*% cbind(FP)-TRP[c(t,(t+1))])))
     b.drawn<-mvrnorm(1, bmean1, b1)
          }else{
     XP<-rbind(Fix[1,], Fix[2,])
     WP<-rbind(Rand[1,], Rand[2,])
     YP<-rbind(Ldep[1], Ldep[2])
     bterm1<- sumtrans(Wstar)
     b1<-pseudoinverse((t(WP)%*% TSig %*%(WP)+bterm1+solve(DDi)))
     b1 <-  make.positive.definite(b1)
     b1 <- matrix(as.real(b1), nrow=RRB)
     Ldephat<-as.matrix(Poly(Rho=phig, XY=as.matrix(Ldep), p=Lag))
     bterm2<-crosssum(XX=Wstar, YY=(Ldephat-Xstar%*%cbind(FP)-cstar))
     bmean1<-t(b1%*%(t(WP)%*%TSig%*%(YP-XP%*% cbind(FP)-TRP[ c(t,(t+1))])+as.matrix(bterm2)))
     b.drawn<-mvrnorm(1, bmean1, b1)
   }
        }
     b.drawn
}

# rrows=rows;Rand=Wi; BSig=Sig; DDi=Di; Fix=Xi; Ldep=ystar; Bccin=Tccin;  Xstar=Xstar; Wstar=Wstar; cstar=cstar; t=min(T) 
BetaTerm <- function (rrows=rows,Rand=Wi, BSig=Sig, DDi=Di, Fix=Xi, Ldep=ystar, Bccin=Tccin,  Xstar=Xstar, Wstar=Wstar, cstar=cstar, p=p, pphig=phig ){
      qq <- matrix(NA, nrow=2, ncol=2)
      qq1 <- c()
      CC<-ncol(DDi)
   if (length(rrows)==1){
         BB1 <- rbind(Rand)%*%DDi%*%cbind(Rand)+BSig[1,1]
          qq<-as.vector(1/BB1)*cbind(Fix)%*%rbind(Fix)
          qq1 <- t(as.matrix(as.vector(1/BB1)*rbind(Fix)* (Ldep-Bccin[1])))
       }else{
           XP<-rbind(Fix[1,], Fix[2,])
         WP<-rbind(Rand[1,], Rand[2,])
         YP<-rbind(Ldep[1], Ldep[2])
           BB1<- WP%*%DDi%*%t(WP)+BSig
          BB11 <- pseudoinverse(BB1)
         if (length(rrows)==2){
        qq <- t(XP)%*%BB11%*%XP
        qq1 <- t(XP)%*%BB11%*%(YP-Bccin[ 1:p])
           }else{
        BB2<-Wstar%*%DDi%*% t(Wstar)+diag(1, (length(Ldep)-p))
          BB12 <- pseudoinverse(BB2)
        Ldephat<-as.matrix(Poly(Rho=pphig, XY=as.matrix(Ldep), p=p))    
        qq<-t(Xstar)%*%BB12%*%Xstar+t(XP)%*%BB11%*%XP
        qq1 <- t(XP)%*%BB11%*%(YP-Bccin[1:p])+(t(Xstar)%*% BB12%*%(Ldephat-cstar))
            }
      }
     return(cbind(qq, qq1))
}

# Function: create diagonal block matrix
bdiag <- function(x){
     if(!is.list(x)) stop("x not a list")
     n <- length(x)
     if(n==0) return(NULL)
     x <- lapply(x, function(y) if(length(y)) as.matrix(y) else
stop("Zero-length component in x"))
     d <- array(unlist(lapply(x, dim)), c(2, n))
     rr <- d[1,]
     cc <- d[2,]
     rsum <- sum(rr)
     csum <- sum(cc)
     out <- array(0, c(rsum, csum))
     ind <- array(0, c(4, n))
     rcum <- cumsum(rr)
     ccum <- cumsum(cc)
     ind[1,-1] <- rcum[-n]
     ind[2,] <- rcum
     ind[3,-1] <- ccum[-n]
     ind[4,] <- ccum
     imat <- array(1:(rsum * csum), c(rsum, csum))
     iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2],
(y[3]+1):y[4]], imat=imat)
     iuse <- as.vector(unlist(iuse))
     out[iuse] <- unlist(x)
     return(out)
} 

# Function to prepare for drawing bivariate normal, return the mean and variance of conditional distibution of each
BiNorm1<-function(mu, y2,Sigma){
  mean1<- mu[1]+Sigma[1,2]*solve(Sigma[2,2])*(y2-mu[2])
  sigma<-Sigma[1,2]*solve(Sigma[2,2])*Sigma[2,1]
  return(c(mean1, sigma))
}

BiNorm2<-function(mu, y1,Sigma){
  mean2<- mu[2]+Sigma[2,1]*solve(Sigma[1,1])*(y1-mu[1])
  sigma2<-Sigma[2,1]*solve(Sigma[1,1])*Sigma[1,2]
  return(c(mean2, sigma2))
}

# function 3: draw from truncated normal from (-inf 0) or its complement
BTnorm<-function(x, mu, sigma){
           upper <- pnorm(0, mu, sigma)
            if(x==0) {
                 u<-runif(1)
                xstar<-qnorm(u*upper)
                 } else{
                 u <- runif(1)
                xstar<-qnorm(u-u*upper+upper)
                 }
         xstar
     }

BTnormlog <- function(x, mu, sigma=1){
            if(x==0) {
                 u <- runif(1)
                xstar <- qnorm(log(u)+pnorm(0, mu, sigma, log=T), mu, sigma, log=T)
                 } else{
                 u <- runif(1)
                xstar <- -qnorm(log(u)+pnorm(0, -mu, sigma, log=T), -mu, sigma, log=T)
                 }
           xstar
        }
#_______________________________________________________________________
# A function creating AR(1) Covariance matrix
Cov.AR1 <- function(rho=rho, col=T, sigma=sigma){
  cov.vector <- sigma
 for (k in 2: col){
   cov.vector[k] <- rho*cov.vector[k-1]
 }
  covariance <- matrix(NA, nrow=col, ncol=col)
  covariance[1,1] <- sigma
  for (i in 1:col){
    for (j in 1:col){
      covariance[i, j] <- cov.vector[abs(i-j)+1]
    }
  }
   covariance
}

#------------------------------------------------------------------------
Cov.AR2 <- function(phig=phig, col=T){
II<-diag(4)
 e<-matrix(c(1,0), nrow=2)
  G<-matrix(c(phig[1],1, phig[2], 0), ncol=2)
  d<-as.vector(pseudoinverse(II-kronecker(G,G))%*%as.vector(e%*%t(e)))
  Sig<-matrix(d, nrow=2)
  var <- Sig[1,1]
  cov.vector <- c(Sig[1,1], Sig[1,2])
 for (k in 3: col){
   cov.vector[k] <- phig[1]*cov.vector[k-1]+phig[2]*cov.vector[k-2]
 }
  covariance <- matrix(NA, nrow=col, ncol=col)
  covariance[1:2, 1:2] <- Sig
  for (i in 1:col){
    for (j in 1:col){
      covariance[i, j] <- cov.vector[abs(i-j)+1]
    }
  }
  return(list(Sig, covariance))
}

#_________________________________________________________________________

Cov.ARp <- function(phig=phig, col=T){
  ord <- length(phig)
  II<-diag(ord^2)
  om <- rep(0, (ord-1))
 e<-matrix(c(1,om), nrow=ord)
  G<-matrix(NA, ncol=ord, nrow=ord)
  G[2:ord, 1:(ord-1)] <- diag(ord-1)
  G[1,] <- phig
  G[-1, ord] <- 0

  d<-as.vector(pseudoinverse(II-kronecker(G,G))%*%as.vector(e%*%t(e)))
  Sig<-matrix(d, nrow=ord)
  var <- Sig[1,1]
  cov.vector <- as.vector(Sig[,1])
 for (k in (ord+1): col){
   cov.vector[k] <- rev(phig) %*% cov.vector[(k-ord): (k-1)]
 }
  covariance <- matrix(NA, nrow=col, ncol=col)
  covariance[1:ord, 1:ord] <- Sig
  for (i in 1:col){
    for (j in 1:col){
      covariance[i, j] <- cov.vector[abs(i-j)+1]
    }
  }
  return(list(Sig, covariance))
}
  
# function: sum the term of X%*% Y for all the observations --- for use of posterior mean (normal) update
crosssum<-function(XX, YY){
     total1 <- rep(0,ncol(XX))
       for (d in 1: nrow(XX)){
            total1<-total1+XX[d,]*YY[d]
       }
     total1
   }


#Fix=Xi;Rand=Wi;T=T;Dep=yi;Ldep=ystari;FP=betan;RP=bbbn;TRP=ccin;unit=n;Lag=p;phig=phig;Xstar=Xstar;Wstar=Wstar;cstar=cstar

CutPointUpdate <- function(obs, latent, oldcut, cut0){
     ntau<-c()
  catnum <- length(oldcut)+2        # number of categories
  if (catnum==3){
    al<-max(max(latent[which(obs==2)]), cut0)
      bl<-min(latent[which(obs==3)])
      ntau<-runif(1, al, bl)
   }
     if (catnum==4){
      a1<-max(max(latent[which(obs==2)]), cut0) 
      b1<-min(min(latent[which(obs==3)]), oldcut[2])
      ntau[1]<-runif(1, a1, b1)
      al<-max(max(latent[which(obs==3)]), ntau[1])
      bl<-min(latent[which(obs==4)])
      ntau[2]<-runif(1, al, bl)
      }
     if (catnum>4){
      a1<-max(max(latent[which(obs==2)]), cut0) 
      b1<-min(min(latent[which(obs==3)]), oldcut[2])
      ntau[1]<-runif(1, a1, b1)
   for (i in 3:(catnum-2) ){
      aaa<-max(max(latent[which(obs==(i))]), ntau[i-2])
      bbb<-min(min(latent[which(obs==(i+1))]), oldcut[i])
      ntau[i-1]<-runif(1, aaa, bbb)
    }
      al<-max(max(latent[which(obs==(catnum-1))]), ntau[catnum-3])
      bl<-min(latent[which(obs==catnum)])
      ntau[catnum-2]<-runif(1, al, bl)
    }
      ntau
 }

Poly<-function(Rho, XY, p){   # XY is matrix, Rho is a vector
  trans<-XY[1,]
  for (tt in (p+1): nrow(XY)){
    trans1<-rep(0, ncol(XY))
  for (w in 1:p){
  trans1<-trans1+XY[tt-w, ]*Rho[w]
   }
  trans<-rbind( trans, rbind(XY[tt,]-trans1))
  }
  rtrans<-trans[-1,]
}


Poly2<-function(Rho, XY, p){   # XY is matrix, Rho is a vector
  trans<-XY[1]
  for (tt in (p+1): length(XY)){
    trans1<-0
  for (w in 1:p){
  trans1<-trans1+XY[tt-w]*Rho[w]
   }
  trans<-rbind( trans, rbind(XY[tt]-trans1))
  }
  rtrans<-trans[-1]
}

#Fix=Xi; Rand=Wi; T=TT; Dep=yi; Ldep=ystari; FP=betan; RP=bbbn; TRP=ccin; unit=n; Lag=p; phig=phig; Xstar=Xstar; Wstar=Wstar; cstar=cstar; Entry=entry; Id1=id1
DataAug<-function(Fix=Xi, Rand=Wi, T=T, Dep=yi, Ldep=ystari, FP=betan, RP=bbbn, TRP=ccin, unit=n, Lag=p, phig=phig, Xstar=Xstar, Wstar=Wstar, cstar=cstar, Entry=entry, Id1=id1){
   ystar <- c()
   rows<-Entry[which(Id1==unit)]
    t<-min(T)
   yl<-length(Dep)
          if (yl==1) {
      mu1 <- Fix %*% cbind(FP)+Rand%*%cbind(RP)+TRP[t]
      sig1 <- sqrt(Sig[1,1])
      ystar[1]=BTnorm(x=Dep[1], mu=mu1, sigma=sig1)
    }else{
          if (yl==2){
      mu12<-c(Fix[1,]%*%cbind(FP)+Rand[1,]%*%cbind(RP)+TRP[t], Fix[2,]%*%cbind(FP)+Rand[2,]%*%cbind(RP)+TRP[t+1])
      x1norm<-BiNorm1(mu=mu12, Sigma=Sig, y2=Ldep[2])
      mu11<-x1norm[1]
      sig<-sqrt(x1norm[2])
       ystar[1]=BTnorm(x=Dep[1], mu=mu11, sigma=sig)
      
      x2norm<-BiNorm2(mu=mu12, Sigma=Sig, y1=ystar[1])
      mu21<-x2norm[1]
      sig2<-sqrt(x2norm[2])
      ystar[2]=BTnorm(x=Dep[2], mu=mu21, sigma=sig2)
       }else{     ### length(rows)>=3
     mu12<-c(Fix[1,]%*%cbind(FP)+Rand[1,]%*%cbind(RP)+TRP[t], Fix[2,]%*%cbind(FP)+Rand[2,]%*%cbind(RP)+TRP[t+1])
      x1norm<-BiNorm1(mu=mu12, Sigma=Sig, y2=Ldep[2])
      mu11<-x1norm[1]
      var1<-1/(x1norm[2])
   
      mu22<-Ldep[3]-phig[1]*Ldep[2]- (Xstar[1,]%*% cbind(FP)+Wstar[1,]%*%cbind(RP)+cstar[1])
   
      sig1<-1/sqrt(var1+phig[2]^2)
      mu1<-sig1^2*(phig[2]*mu22+var1*mu11)
   
      ystar[1]=BTnorm(x=Dep[1], mu=mu1, sigma=sig1)

    ################## y_{i2} ###################################################
   
         x2norm<-BiNorm2(mu=mu12, Sigma=Sig, y1=ystar[1])
         mu12<-x2norm[1]
         var2 <-1/(x2norm[2])
   
         mu23<- Ldep[3] -phig[2] * ystar[1]-(Xstar[1,]%*% cbind(FP)+Wstar[1,]%*%cbind(RP)+cstar[1])
         if (yl==3){
            sig2<-1/sqrt(phig[1]^2+var2)
         mu2<-sig2^2*(phig[1]*mu23 +var2*mu12)
   
         ystar[2]<-BTnorm(x=Dep[2], mu=mu2, sigma=sig2)
          }else{
          mu24<-Ldep[4]- phig[1] * Ldep[3]-(Xstar[2,]%*% cbind(FP)+Wstar[2,]%*%cbind(RP)+cstar[2] )
   
         sig2<-1/sqrt(phig[1]^2+phig[2]^2+var2)
         mu2<-sig2^2*(phig[1]*mu23 +phig[2]*mu24+var2*mu12)
   
         ystar[2]<-BTnorm(x=Dep[2], mu=mu2, sigma=sig2)
        }
}
}
     
################  y_{it} ##################################################
      if (yl>=3){
      for (i in 3:yl) {
          if (i < (yl-1)){     #### ith unit has at least 5 observations
           mu4 <-Xstar[i-Lag,]%*% cbind(FP)+Wstar[i-Lag,]%*%cbind(RP)+cstar[i-Lag] + phig[1] * ystar[i-1]+phig[2] * ystar[i-2]
           
            mu5 <- Ldep[i+1] -phig[2] * ystar[i-1]-(Xstar[i-Lag+1,]%*% cbind(FP)+Wstar[i-Lag+1,]%*%cbind(RP)+cstar[i-Lag+1]) 
           
            mu6 <-Ldep[i+2]- phig[1] * Ldep[i+1]-(Xstar[i-Lag+2,]%*% cbind(FP)+Wstar[i-Lag+2,]%*%cbind(RP)+cstar[i-Lag+2] )
           
           vmu<-1/sqrt(phig[1]^2+phig[2]^2+1)
           tmu<-vmu^2*(phig[1]*mu5 +phig[2]*mu6+mu4)

           ystar[i]=BTnorm(x=Dep[i], mu=tmu, sigma=vmu)
         }else{
           if (i == (yl-1)){
            mu4 <-Xstar[i-Lag,]%*% cbind(FP)+Wstar[i-Lag,]%*%cbind(RP)+cstar[i-Lag] + phig[1] * ystar[i-1]+phig[2] * ystar[i-2]
           
            mu5 <- Ldep[i+1] -phig[2] * ystar[i-1]-(Xstar[i-Lag+1,]%*% cbind(FP)+Wstar[i-Lag+1,]%*%cbind(RP)+cstar[i-Lag+1]) 
           
              vmu<-1/sqrt(1+phig[1]^2)
             tmu<-vmu^2*(phig[1]*mu5 +mu4)
            
              ystar[i]<-BTnorm(x=Dep[i], mu=tmu, sigma=vmu)
            
           }else{
               mu4 <-Xstar[i-Lag,]%*% cbind(FP)+Wstar[i-Lag,]%*%cbind(RP)+cstar[i-Lag] + phig[1] * ystar[i-1]+phig[2] * ystar[i-2]
           
           ystar[i]=BTnorm(x=Dep[i], mu=mu4, sigma=1)
           }
        }
        }
}
    ystar
}

############## Data Augmentation by using polynomial transformation ##############


DataAug2<-function(Fix=Xi, Rand=Wi, T=T, Dep=yi, Ldep=ystari, FP=betan, RP=bbbn, TRP=rant, unit=n, Lag=p, phig=phig, Sig=Sig, Entry=entry, Id1=id1, categ=categ, tau, Y=Y){
   ystar <- c()
   rows<-Entry[which(Id1==unit)]
   yl<-length(Dep)
   if (yl==1) {
      mu1 <- Fix %*% cbind(FP)+Rand%*%cbind(RP)+TRP
      sig1 <- sqrt(Sig[1,1])
      if (categ==2){
      ystar[1] <- BTnorm(x=Dep[1], mu=mu1, sigma=sig1)
}else{
      ystar[1] <- OrdinalAugment(obs=Dep[1], mu=mu1, sigma=sig1, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
}
    }else{
          if (yl==2){
      mu12<-c(Fix[1,]%*%cbind(FP)+Rand[1,]%*%cbind(RP)+TRP[1], Fix[2,]%*%cbind(FP)+Rand[2,]%*%cbind(RP)+TRP[2])
      x1norm<-BiNorm1(mu=mu12, Sigma=Sig, y2=Ldep[2])
      mu11<-x1norm[1]
      sig<-sqrt(x1norm[2])
       ystar[1]=BTnorm(x=Dep[1], mu=mu11, sigma=sig)
      
      x2norm<-BiNorm2(mu=mu12, Sigma=Sig, y1=ystar[1])
      mu21<-x2norm[1]
      sig2<-sqrt(x2norm[2])
      if (categ==2){
      ystar[2] <- BTnorm(x=Dep[2], mu=mu21, sigma=sig2)
    }else{
       ystar[2] <- OrdinalAugment(obs=Dep[2], mu=mu21, sigma=sig2, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
    }
       }else{     ### length(rows)>=3
       Xstar <- Poly(Rho=phig, XY=Fix, p=Lag)
       Xstar <- rbind(Fix[1:Lag,], Xstar)
       Wstar <- Poly(Rho=phig, XY=Rand, p=Lag)
       Wstar <- rbind(Rand[1:Lag,], Wstar)
       Cstar <- Poly2(Rho=phig, XY=TRP, p=Lag)
       Cstar <- c(TRP[1:Lag], Cstar)
      mu12<-c(Fix[1,]%*%cbind(FP)+Rand[1,]%*%cbind(RP)+TRP[1], Fix[2,]%*%cbind(FP)+Rand[2,]%*%cbind(RP)+TRP[2])
      x1norm<-BiNorm1(mu=mu12, Sigma=Sig, y2=Ldep[2])
      mu11<-x1norm[1]
      var1<-1/(x1norm[2])
   
      mu22<-Ldep[3]-phig[1]*Ldep[2]- (Xstar[3,]%*% cbind(FP)+Wstar[3,]%*%cbind(RP)+Cstar[3])
   
      sig1<-1/sqrt(var1+phig[2]^2)
      mu1<-sig1^2*(phig[2]*mu22+var1*mu11)
     if (categ==2){
      ystar[1] <- BTnorm(x=Dep[1], mu=mu1, sigma=sig1)
}else{
      ystar[1] <- OrdinalAugment(obs=Dep[1], mu=mu1, sigma=sig1, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
}
     

    ################## y_{i2} ###################################################
   
         x2norm<-BiNorm2(mu=mu12, Sigma=Sig, y1=ystar[1])
         mu12<-x2norm[1]
         var2 <-1/(x2norm[2])
   
         mu23<- Ldep[3] -phig[2] * ystar[1]-(Xstar[3,]%*% cbind(FP)+Wstar[3,]%*%cbind(RP)+Cstar[3])
         if (yl==3){
            sig2<-1/sqrt(phig[1]^2+var2)
         mu2<-sig2^2*(phig[1]*mu23 +var2*mu12)
   
         ystar[2]<-BTnorm(x=Dep[2], mu=mu2, sigma=sig2)
          }else{
          mu24<-Ldep[4]- phig[1] * Ldep[3]-(Xstar[4,]%*% cbind(FP)+Wstar[4,]%*%cbind(RP)+Cstar[4] )
   
         sig2<-1/sqrt(phig[1]^2+phig[2]^2+var2)
         mu2<-sig2^2*(phig[1]*mu23 +phig[2]*mu24+var2*mu12)
    if (categ==2){
      ystar[2] <- BTnorm(x=Dep[2], mu=mu21, sigma=sig2)
    }else{
       ystar[2] <- OrdinalAugment(obs=Dep[2], mu=mu21, sigma=sig2, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
    }
        }
}
}

     
################  y_{it} ##################################################
      if (yl>=3){
       Xstar <- Poly(Rho=phig, XY=Fix, p=Lag)
       Xstar <- rbind(Fix[1:Lag,], Xstar)
       Wstar <- Poly(Rho=phig, XY=Rand, p=Lag)
       Wstar <- rbind(Rand[1:Lag,], Wstar)
       Cstar <- Poly2(Rho=phig, XY=TRP, p=Lag)
       Cstar <- c(TRP[1:Lag], Cstar)
      for (i in 3:yl) {
          if (i < (yl-1)){     #### ith unit has at least 5 observations
           mu4 <-Xstar[i,]%*% cbind(FP)+Wstar[i,]%*%cbind(RP)+Cstar[i] + phig[1] * ystar[i-1]+phig[2] * ystar[i-2]
           
            mu5 <- Ldep[i+1] -phig[2] * ystar[i-1]-(Xstar[i+1,]%*% cbind(FP)+Wstar[i+1,]%*%cbind(RP)+Cstar[i+1]) 
           
            mu6 <-Ldep[i+2]- phig[1] * Ldep[i+1]-(Xstar[i+2,]%*% cbind(FP)+Wstar[i+2,]%*%cbind(RP)+Cstar[i+2] )
           
           vmu<-1/sqrt(phig[1]^2+phig[2]^2+1)
           tmu<-vmu^2*(phig[1]*mu5 +phig[2]*mu6+mu4)
             if (categ==2){
     ystar[i]=BTnorm(x=Dep[i], mu=tmu, sigma=vmu)
}else{
      ystar[i] <- OrdinalAugment(obs=Dep[i], mu=tmu, sigma=vmu, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
}
          
         }else{
           if (i == (yl-1)){
            mu4 <-Xstar[i,]%*% cbind(FP)+Wstar[i,]%*%cbind(RP)+Cstar[i] + phig[1] * ystar[i-1]+phig[2] * ystar[i-2]
           
            mu5 <- Ldep[i+1] -phig[2] * ystar[i-1]-(Xstar[i+1,]%*% cbind(FP)+Wstar[i+1,]%*%cbind(RP)+Cstar[i+1]) 
           
              vmu<-1/sqrt(1+phig[1]^2)
             tmu<-vmu^2*(phig[1]*mu5 +mu4)
            
                 if (categ==2){
     ystar[i]=BTnorm(x=Dep[i], mu=tmu, sigma=vmu)
}else{
      ystar[i] <- OrdinalAugment(obs=Dep[i], mu=tmu, sigma=vmu, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
}
            
           }else{
               mu4 <-Xstar[i,]%*% cbind(FP)+Wstar[i,]%*%cbind(RP)+Cstar[i] + phig[1] * ystar[i-1]+phig[2] * ystar[i-2]
           
                if (categ==2){
     ystar[i]=BTnorm(x=Dep[i], mu=mu4, sigma=1)
}else{
      ystar[i] <- OrdinalAugment(obs=Dep[i], mu=mu4, sigma=1, maxcat=max(Y), mincat=min(Y), cutpoint=tau)
}
           }
        }
        }
 }
    ystar
}
######## Data Augmentation by using Geweke method ###########################
DataAug3 <- function(yold, ymean, yi, Cov, categ=categ, tau=tau[g-1,], Y=Y){
  ob <- length(yi)
  ystar <- c()
  for (w in 1:ob){
    u1 <- ymean[w]
    u2 <- ymean[-w]
    var1 <- Cov[w,w]
    if (ob==2){
    var2 <- Cov[-w,][-w]
}else{
    var2 <- Cov[-w,][,-w]
}
    cov12 <- rbind(Cov[w,][-w])
    cov21 <- cbind(Cov[,w][-w])
    if (w==1){
    x2 <-yold[(w+1):ob]
     }else{
    if (w==ob){
    x2 <- ystar[1:(w-1)]
  }else{
    x2 <- c(ystar[1:(w-1)], yold[(w+1):ob])
    }
}
    mean <- u1+cov12%*%solve(var2)%*%(x2-u2)
    var <- var1-cov12%*%solve(var2)%*%cov21
    if (categ==2){
    ystar[w] <- BTnormlog(x=yi[w], mu=mean, sigma=sqrt(var))
}else{
    ystar[w] <- OrdinalAugment(obs=yi[w], mu=mean, sigma=sqrt(var), maxcat=max(Y), mincat=min(Y), cutpoint=tau)
}
}
  return(ystar)
}

#************************* Function to arrange the data **********************#
# DataArrange --- a function to put the covariates at different levels to be  #
#                  at the single level, ready for analysis by the main func.  #
#                                                                             #
# Usage --- id1: index of unit                                                #
#           id2: index of time                                                #
#             y: response variable                                            #
#           Fixed: Covariates of fixed effects                                #
#           UnitRandom: Covariates of unit-specific random effects            #
#           TimeRandom: Covariates of time-specific random effects            #
#           UnitPred: unit-level predictors                                   #
#           TimePred: time-level predictors                                   #
#           rescale: standardize the covariates? 0--FALSE, 1--TRUE            #
#           cont: continuous responses?   0---FALSE, 1---TRUE                 #
#           if there is no covariates, put "NULL"
#                                                                             #
# Data of Creation: 1/1/2009                                                  #
#                                                                             #
#*****************************************************************************#
# This function is for the model specification with group level predictors at both unit and time group level
GroupPred2 <- function(X, index){
  b<-ncol(X)
  J <- length(unique(na.omit(index)))
  XPred<-matrix(0, nrow=J, ncol=b)
  for (i in 1:b){
    if(length(which(is.na(X[,i])))>0){
      average <- tapply(na.omit(X[,i]), index[-which(is.na(X[,i]))], mean)
      XPred[, i] <- average
    }else{
      average <- tapply(X[,i], index, mean)
      XPred[, i] <- average
    }
}
}

DataArrange <- function(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred, TimePred=TimePred){
  # Make sure the design matrices are matrix formats
   if(!is.matrix(Fixed)){Fixed <- as.matrix(Fixed)}
   if(!is.matrix(UnitRandom)){UnitRandom <- as.matrix(UnitRandom)}
   if(!is.matrix(TimeRandom)){TimeRandom <- as.matrix(TimeRandom)}
   if(!is.matrix(UnitPred)){UnitPred <- as.matrix(UnitPred)}
   if(!is.matrix(TimePred)){TimePred <- as.matrix(TimePred)}

   J<-length(unique(na.omit(id1)))  # how many units, in case NAs
   len <- 1:length(na.omit(id1))    # how many observations in the dataset

#________________________ Create group-level predictors____________________#

 # create a design matrix  at the group level of units
   UnitP <-GroupPred(X=as.matrix(UnitPred), index=id1)
   A <- as.matrix(UnitP)      
 # create a design matrix  at the group level of units
   UnitQ <-GroupPred(X=as.matrix(TimePred), index=id2) 
   B <- as.matrix(UnitQ)   

#_______________________ Interacting terms in the individual level _______#

# interacting terms of unit-level variables and individual variables--- they are with fixed effects in the reduced model
   RanFix <- random(Var=UnitRandom, Pred=A, index=id1)

 # interacting terms of time-level variables and individual variables--- they are with fixed effects in the reduced model
   TimFix <- randomt2(Var=TimeRandom, Pred=B, index1=id1, index2=id2, entry=len)

#_______________________ Get the dimensions _______________________________#

   RD <- ncol(RanFix)
   TD <- ncol(TimFix)
   FD <- ncol(Fixed)
   UD <- ncol(UnitRandom)
   UT <- ncol(TimeRandom)

#__________________________ Creat the dataset ______________________________#

# conbined the data with two indices
   RawData<-cbind(id1, id2, y, Fixed, RanFix, TimFix, UnitRandom, TimeRandom) 
   Data2 <- remove.NA(RawData)

 # get rid of entries with missing data---no drop-out: THIS NEEDS FURTHER THINKING---HOW TO DEAL WITH MISSING DATA

   DD <- ncol(Data2)
#_______________________ Recale the data and get right covariates ___________#
   newid1 <- Data2[,1]
   newid2 <- Data2[,2]
   Y<-Data2[,3]
   X <- as.matrix(cbind(Data2[, c(4: (3+FD+TD+RD))]))
   W <- as.matrix(Data2[,  c((4 +FD+TD+RD):(3+FD+TD+RD+UD))])
   S <- as.matrix(Data2[, (4+FD+TD+RD+UD):(3+FD+TD+RD+UD+UT)])

   return(list(newid1, newid2,Y, X, W, S))
}

# This is the function only have unit-level predictors. If time-specific covariates is only an intercept, the intercept is dropped from X, only a time-varing intercept will be estimated. If universe intercept is needed, it has to be specified in Fixed
DataArrange.UP <- function(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred){
  # Make sure the design matrices are matrix formats
   if(!is.matrix(Fixed)){Fixed <- as.matrix(Fixed)}
   if(!is.matrix(UnitRandom)){UnitRandom <- as.matrix(UnitRandom)}
   if(!is.matrix(TimeRandom)){TimeRandom <- as.matrix(TimeRandom)}
   if(!is.matrix(UnitPred)){UnitPred <- as.matrix(UnitPred)}

   J<-length(unique(na.omit(id1)))  # how many units, in case NAs
   len <- 1:length(na.omit(id1))    # how many observations in the dataset

#________________________ Create group-level predictors______________#

 # create a design matrix  at the group level of units
   UnitP <-GroupPred(X=as.matrix(UnitPred), index=id1)
   A <- as.matrix(UnitP)      
#_______________________ Interacting terms in the individual level _____#
# interacting terms of unit-level variables and individual variables--they are with fixed effects in the reduced model
   RanFix <- random(Var=UnitRandom, Pred=A, index=id1)
#_______________________ Get the dimensions _________________________#

   RD <- ncol(RanFix)
   FD <- ncol(Fixed)
   UD <- ncol(UnitRandom)
   UT <- ncol(TimeRandom)

#__________________________ Creat the dataset _______________________
# conbined the data with two indices
  if (ncol(TimeRandom)==1&length(unique(TimeRandom))==1){
    RawData<-cbind(id1, id2, y, Fixed,  RanFix, UnitRandom)
  }else{
   RawData<-cbind(id1, id2, y, Fixed, RanFix, UnitRandom, TimeRandom)
}
   Data2 <- remove.NA(RawData)

 # get rid of entries with missing data---no drop-out: THIS NEEDS FURTHER THINKING---HOW TO DEAL WITH MISSING DATA 
   DD <- ncol(Data2)
#_______________________ Recale the data and get right covariates ___

   newid1 <- Data2[,1]
   newid2 <- Data2[,2]
   Y<-Data2[,3]
   X <- as.matrix(cbind(Data2[, c(4: (3+FD+RD))]))
   W <- as.matrix(Data2[, c((4+FD+RD):(3+FD+RD+UD))])
   if (ncol(TimeRandom)==1&length(unique(TimeRandom))==1){
     S <- as.matrix(rep(1, nrow(X)))
   }else{
     S <- as.matrix(Data2[, c((4+FD+RD+UD):(3+FD+RD+UD+UT))])
   }
   return(list(newid1, newid2,Y, X, W, S))
}

# This is the function only have unit-level predictors. If unit-specific covariates is only an intercept, the intercept is dropped from X, only a unit-varing intercept will be estimated. If universe intercept is needed, it has to be specified in Fixed
DataArrange.TP <- function(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom,  TimePred=TimePred){
  # Make sure the design matrices are matrix formats
   if(!is.matrix(Fixed)){Fixed <- as.matrix(Fixed)}
   if(!is.matrix(UnitRandom)){UnitRandom <- as.matrix(UnitRandom)}
   if(!is.matrix(TimeRandom)){TimeRandom <- as.matrix(TimeRandom)}
   if(!is.matrix(TimePred)){TimePred <- as.matrix(TimePred)}

   J<-length(unique(na.omit(id1)))  # how many units, in case NAs
   len <- 1:length(na.omit(id1))    # how many observations in the dataset

#________________________ Create group-level predictors____________________#
 # create a design matrix  at the group level of units
   UnitQ <-GroupPred(X=as.matrix(TimePred), index=id2) 
   B <- as.matrix(UnitQ)   

#_______________________ Interacting terms in the individual level _______#
 # interacting terms of time-level variables and individual variables--- they are with fixed effects in the reduced model
   TimFix <- randomt2(Var=TimeRandom, Pred=B, index1=id1, index2=id2, entry=len)

#_______________________ Get the dimensions _______________________________#

   TD <- ncol(TimFix)
   FD <- ncol(Fixed)
   UD <- ncol(UnitRandom)
   UT <- ncol(TimeRandom)

#__________________________ Creat the dataset ______________________________#

   if (ncol(UnitRandom)==1&length(unique(UnitRandom))==1){
    RawData<-cbind(id1, id2, y, Fixed,   TimeRandom)
  }else{
# conbined the data with two indices
   RawData<-cbind(id1, id2, y, Fixed, UnitRandom, TimeRandom, TimFix)
 }
   Data2 <- remove.NA(RawData)

 # get rid of entries with missing data---no drop-out: THIS NEEDS FURTHER THINKING---HOW TO DEAL WITH MISSING DATA
   DD <- ncol(Data2)
   DDT <- ncol(Data)
#_______________________ Recale the data and get right covariates ___________#
   newid1 <- Data2[,1]
   newid2 <- Data2[,2]
   Y<-Data2[,3]
   X <- as.matrix(cbind(Data2[, -c(c(1:3), c((DDT-UT+1):DDT))]))
   if (ncol(UnitRandom)==1&length(unique(UnitRandom))==1){
     W <- as.matrix(rep(1, nrow(X)))
}else{
   W <-  as.matrix(Data2[, (4+FD):(3+FD+UD)])
}
   S <- as.matrix(Data2[,  c((DDT-UT+1):DDT)])

   return(list(newid1, newid2,Y, X, W, S))
}

# Function: there is no group level predictors, if only varing-intercepts are included, the universal intercept should not be included.
DataArrange.NP <- function(id1=Unit, id2=Time, y=y, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom){
    if(!is.matrix(Fixed)){Fixed <- as.matrix(Fixed)}
   if(!is.matrix(UnitRandom)){UnitRandom <- as.matrix(UnitRandom)}
   if(!is.matrix(TimeRandom)){TimeRandom <- as.matrix(TimeRandom)}

  # Make sure the design matrices are matrix formats
   J<-length(unique(na.omit(id1)))  # how many units, in case NAs
   len <- 1:length(na.omit(id1))    # how many obsedrvations in the dataset

   FD <- ncol(Fixed)
   UD <- ncol(UnitRandom)
   UT <- ncol(TimeRandom)
#__________________________ Creat the dataset ______________________________#

# conbined the data with two indices
   if ((UD==1&length(unique(UnitRandom))==1)&(UT==1&length(unique(TimeRandom))==1)){
     RawData<-cbind(id1, id2, y, Fixed)
     Data <- remove.NA(RawData)
     newid1 <- Data[,1]
     newid2 <- Data[,2]
      Y<-Data[,3]
      X <- as.matrix(Data[, -c(1:3)])
      W <- as.matrix(rep(1,nrow(X)))
      S <- as.matrix(rep(1,nrow(X)))
   }else{
    if ((UD==1&length(unique(UnitRandom))==1)&(UT!=1|length(unique(TimeRandom))!=1)){
       RawData<-cbind(id1, id2, y, Fixed,  TimeRandom)
       Data <- remove.NA(RawData)
       RR <- ncol(Data)
        newid1 <- Data[,1]
        newid2 <- Data[,2]
        Y<-Data[,3]
        X <- as.matrix(Data[, -c(1:3)])
        W <- as.matrix(rep(1,nrow(X)))
        S <- as.matrix(Data[, (RR-UT+1):RR])
     }else{
        if ((UD!=1|length(unique(UnitRandom))!=1)&(UT==1&length(unique(TimeRandom))==1)){
          RawData<-cbind(id1, id2, y, Fixed, UnitRandom)

       Data <- remove.NA(RawData)
       RR <- ncol(Data)
        newid1 <- Data[,1]
        newid2 <- Data[,2]
        Y<-Data[,3]
        X <- as.matrix(Data[, -c(1:3)])
        W <- as.matrix(Data[, (RR-UD+1):RR])
        S <- as.matrix(rep(1,nrow(X)))
        }else{
           RawData<-cbind(id1, id2, y, Fixed, UnitRandom, TimeRandom)
        Data <- remove.NA(RawData)
       RR <- ncol(Data)
        newid1 <- Data[,1]
        newid2 <- Data[,2]
        Y<-Data[,3]
        X <- as.matrix(Data[, -c(1:3)])
        W <- as.matrix(Data[, (RR-UD -UT+1):(RR-UT)])
        S <- as.matrix(Data[, (RR-UT+1):RR])
         }
      }
}

   return(list(newid1, newid2,Y, X, W, S))
}
RegressionData <- function(yob=y, Unit, Time, Fixed, UnitRandom, TimeRandom, UnitPred, TimePred){
    if(!is.numeric(UnitPred)& !is.numeric(TimePred)){
      NewData <- DataArrange.NP(id1=Unit, id2=Time, y=yob, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom)
    } else{
        if (!is.numeric(UnitPred) &is.numeric(TimePred)){
           NewData <- DataArrange.TP(id1=Unit, id2=Time, y=yob, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, TimePred=TimePred)
        } else{
           if (is.numeric(UnitPred) &!is.numeric(TimePred)){
              NewData <- DataArrange.UP(id1=Unit, id2=Time, y=yob, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred)
             }else{
               NewData <- DataArrange(id1=Unit, id2=Time, y=yob, Fixed=Fixed, UnitRandom=UnitRandom, TimeRandom=TimeRandom, UnitPred=UnitPred, TimePred=TimePred)
             }
        }
   }
    return(NewData)
}
DataArrange.ICL <- function(Unit=Unit, UnitPred=UnitPred, UnitRandom=UnitRandom, Fixed=Fixed){
    id<-Index(Unit)                              # create indices of units
    J<-length(unique(na.omit(id)))               # how many units in the dataset
    UnitP <-GroupPred(X=as.matrix(UnitPred), index=id) 
    A <- as.matrix( UnitP)    # create a design matrix  at the group level of units
   RanFix <- random(Var=UnitRandom, Pred=A, index=id)
   RD <- ncol(RanFix)
   FD <- ncol(Fixed)
   UD <- ncol(UnitRandom)
   RawData<-cbind(id, y, Fixed, UnitRandom) # conbined the data with two indices
   Data <- remove.NA(RawData)
   Data2 <- cbind(Data, RanFix)  
   DD <- ncol(Data2)
   Y<-Data2[,2]
   X <- as.matrix(cbind(Data2[, c(3: (2+FD),(UD+3+FD):DD)]))
   W <- as.matrix(Data2[,  c((3 +FD):(UD+2+FD))])
   newdi1 <- Data2[,1]
   return(list(newid1,Y, X, W))
  }



random<-function(Var, Pred, index){
   J<-length(unique(na.omit(index)))
     RanFix <- matrix(0, nrow=1, ncol=ncol(Var)*ncol(Pred))
  for (j in 1:J){
    W <- as.matrix(Var[which(index==j),])
    a <- kronecker(Diagonal(ncol(Var)), t(as.matrix(Pred[j,])))
    RanFix<-rbind(RanFix, as.matrix(W%*%a))
   }
   TTT<-RanFix[-1,]
}


# with indices without NA
randomt2 <- function(Var, Pred, index1, index2, entry){
   index <- na.omit(index1)
    index2 <- na.omit(index2)
   Var <- Var
   Pred <- na.omit(Pred)
   J1<-length(unique(na.omit(index)))
    WWW <- matrix(NA, nrow=length(entry), ncol=ncol(Var)*ncol(Pred))
   for (j in 1:J1){
     rows <- entry[which(index==j)]
     TT <- index2[rows]
     W <- as.matrix(Var[rows,])
     PW <- Pred[TT,]
     Ee <- length(rows)
     Tt <- matrix(NA, ncol=ncol(Var)*ncol(Pred), nrow=Ee)
     for (s in 1:Ee){
       Tt[s,] <- as.vector(t(cbind(W[s,])%*%rbind(PW[s,])))
     }
     WWW[rows,] <- Tt
       
   }
   WWW
}
  


# Function: calculate group level mean
GroupPred<-function(X, index){
  b<-ncol(X)
  J <- length(unique(na.omit(index)))
  XPred<-matrix(0, nrow=J, ncol=b)
  for (i in 1:b){
    Y<-X[,i]
    a <- length(Y)
    Z <- is.na(X[,i])==TRUE
    m <- which(Z==TRUE)
    s <- length(Z[which(Z==FALSE)])
    if (s==a){
    XP <- tapply(Y, index, mean)
    XPred[,i]<-XP
  }else{
    XP <- tapply(X[-m, i], index[-m], mean)
     XPred[,i]<-XP
    }
     XGroup<-XPred
}
}


# function 1: to get  index variable, in which "Names" is the grouping creterion
Index<-function(Names){
uniq <- unique(na.omit(Names))                # get units
J <- length(uniq)                          # how many units 
index <- rep (NA, J)                      # how many observations for each unit
for (i in 1:J) index[Names==uniq[i]] <- i
index
}


myrwish <- function (v=2*FB,S=Dinv ) 
{
    if (!is.matrix(S)) 
        S <- matrix(S)
    if (nrow(S) != ncol(S)) {
        stop(message = "S not square in rwish().\n")
    }
    if (v < nrow(S)) {
        stop(message = "v is less than the dimension of S in rwish().\n")
    }
    p <- nrow(S)
    CC <- sechol(S)
    Z <- matrix(0, p, p)
    diag(Z) <- sqrt(rchisq(p, v:(v - p + 1)))
    if (p > 1) {
        pseq <- 1:(p - 1)
        Z[rep(p * pseq, pseq) + unlist(lapply(pseq, seq))] <- rnorm(p * 
            (p - 1)/2)
    }
    return (crossprod(Z %*% CC))
}
myriwish <- function(v, S)
  {
    return(pseudoinverse(myrwish(v, pseudoinverse(S))))
}

my.positive.making <- function (m, tol) 
{
    if (!is.matrix(m)) 
        m <- as.matrix(m)
    d <- dim(m)[1]
    if (dim(m)[2] != d) 
        stop("Input matrix is not square!")
    es <- eigen(m)
    esv <- as.real(es$values)
    if (missing(tol)) 
        tol <- d * max(abs(esv)) * .Machine$double.eps
    delta <- 2 * tol
    tau <- pmax(0, delta - esv)
    dm <- es$vectors %*% diag(tau, d) %*% t(es$vectors)
    return(m + dm)
}

# obs=yi[j]; mu=mean.ystar[j]; sigma=sqrt(nK); maxcat=max(Y); mincat=min(Y), cutpoint=tau[g-1,]
OrdinalAugment <- function(obs, maxcat, mincat, cutpoint, mu, sigma){
  if(length(cutpoint)!=(maxcat-2)){
    warning("numbers of categories and cutpoints do not match")
        return(NULL)
    }
  ystar <- c()
 if (obs==mincat){
     ystar<-rtnorm(1, mu, sigma,  upper=0)
   }else{
     if  (obs==(mincat+1)){
         ystar<-rtnorm(1, mu, sigma, lower=0,  upper=cutpoint[(obs-1)])
       }else{
     if (obs==maxcat){
       ystar<-rtnorm(1, mu, sigma, lower=cutpoint[length(cutpoint)])
     }else{
       ystar <- rtnorm(1, mu, sigma, lower=cutpoint[(obs-2)], upper=cutpoint[(obs-1)])
     }
   }
   }
  ystar
}
  

rnbnd <- function(N, mn=0, sd=1, lo=-Inf, hi=Inf, tol=3) {
  if (hi-mn < tol*sd || mn-lo < tol*sd) # Alan Zaslavsky's method
    return(mn + sd * qnorm(runif(N, pnorm((lo-mn)/sd), pnorm((hi-mn)/sd))))
   x <- rep(NA, N) # Rejection method
  bnd <- function(z, lo, hi) ifelse(z<lo | z>hi, NA, z)
  while (any(q <- is.na(x))) x[q] <- bnd(rnorm(sum(q), mn, sd), lo, hi)
  x
}

OrdinalAugment2 <- function(obs, maxcat, mincat, cut0, cutpoint, mu, sigma){
  if(length(cutpoint)!=(maxcat-2)){
    warning("numbers of categories and cutpoints do not match")
        return(NULL)
    }
  ystar <- c()
 if (obs==mincat){
     ystar<-rnbnd(1, mn=mu, sd=sigma,   hi=cut0)
   }else{
     if  (obs==(mincat+1)){
         ystar<-rnbnd(1, mn=mu, sd=sigma, lo=cut0,  hi=cutpoint[(obs-1)])
       }else{
     if (obs==maxcat){
       ystar<-rnbnd(1, mn=mu, sd=sigma, lo=cutpoint[length(cutpoint)])
     }else{
       ystar <- rnbnd(1, mn=mu, sd=sigma, lo=cutpoint[(obs-2)], hi=cutpoint[(obs-1)])
     }
   }
   }
  ystar
}

# function: outer product of matrix
 outerself<-function(x) outer(x,x)

#### Good function of tapply
# cty.mns = tapply(y,index,mean) 
# cty.vars = tapply(y,index,var) 
# tapply(X, INDEX, FUN = NULL, ..., simplify = TRUE)
# T=T; Ldep=ystar; Fix=Xi; Rand=Wi; Pbeta=betan; Pbbn=bbbn; ccin=ccin; p=p
PhiTerm <- function(T=T, Ldep=ystar, Fix=Xi, Rand=Wi, Pbeta=betan, Pbbn=bbbn, ccin=ccin, p=p){
     t <- min(T)
     G<-cbind(rep(1, (max(T)-t+1)))
       markystar<-trans2(YY=Ldep, XX=as.matrix(cbind(Fix, Rand)), BB=c(Pbeta,Pbbn))-ccin[ t:max(T)]
       Ymtotal<-0
          Ymtotal<-matrix(0, nrow=2, ncol=2)
        for(q in 3:length(markystar)){
        Ymtotal<-Ymtotal+rbind(markystar[q-1], markystar[q-2])%*%cbind(markystar[q-1], markystar[q-2])
       }
         Ymtotal2<-0
       for(p in 3:length(markystar)){
       Ymtotal2<-Ymtotal2+rbind(markystar[p-1], markystar[p-2])*markystar[p]
       }
       if (kappa(Ymtotal)>10^3){
        nHatPhi<-pseudoinverse(Ymtotal+0.1)
     }else{
        nHatPhi<-pseudoinverse(Ymtotal)
      }
        nhatphi<-nHatPhi%*%Ymtotal2
        return(cbind(nHatPhi, nhatphi))
}

#Entry=entry; unit=id1; NN=GU; yystar=upystar; PX=X; PW=W; Ti=id2; CC=ccin; newbeta=beta[g,]
PhiTerm1 <- function(Entry=entry,  unit=id1, NN=GU, yystar,newb=upbbn, PX=X, PW=W,  newbeta=beta[g,]){
      Ymtotal<-0
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- as.matrix(PX[rows,])
      upW <- as.matrix(PW[rows,])
      marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])
         
      if (length(rows)>1){
         sumt2<-crossprod(marki)-marki[length(marki)]^2
         that<-0
       for (i in 2: length(marki)){
           that<-that+marki[i-1]*marki[i]
       }
         Ymtotal <- Ymtotal + sumt2
         Ymtotal2 <-Ymtotal2+ that
    }
    }
    return (c(Ymtotal, Ymtotal2))
}

#Entry=entry; unit=id1; NN=GU; yystar=upystar; PX=X; PW=W; Ti=id2; CC=ccin; newbeta=beta[g,]
PhiTerm2 <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W, Ti=id2, CC=ccin, newbeta=beta[g,]){
      Ymtotal<-matrix(0, nrow=2, ncol=2)
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- PX[rows,]
      upW <- PW[rows,]
      TT<-Ti[rows]
      upccin <- CC[TT]
      marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])-upccin
         
      if (length(rows)>2){
      for(q in 3:length(marki)){
        Ymtotal<-Ymtotal+rbind(marki[q-1], marki[q-2])%*%cbind(marki[q-1], marki[q-2])
       Ymtotal2<-Ymtotal2+rbind(marki[q-1], marki[q-2])*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}

PhiTerm3 <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W, Ti=id2, CC=ccin, newbeta=beta[g,], lag=lag){
      Ymtotal<-matrix(0, nrow=2, ncol=2)
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- PX[rows,]
      upW <- PW[rows,]
      TT<-Ti[rows]
      upccin <- CC[TT]
      marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])-upccin
        if (length(rows)>lag){
      for(q in (lag+1):length(marki)){
        trans <- c()
        for (w in 1:lag){
          trans[w] <-marki[q-w]
        }
        Ymtotal<-Ymtotal+cbind(trans)%*%rbind(trans)
       Ymtotal2<-Ymtotal2+cbind(trans)*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}

PhiTerm0 <- function(Entry=entry,  unit=id1, NN=GU, yystar,newb=upbbn, PX=X, PW=W, PS=S, Ti=id2, CC=ccin, newbeta=beta[g,]){
      Ymtotal<-0
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- as.matrix(PX[rows,])
      upW <- as.matrix(PW[rows,])
      upS <- as.matrix(PS[rows,])
      TT<-Ti[rows]
      upccin <- as.matrix(CC[TT,])
        JP <- length(rows)
        rant <- c()
      for (d in 1: JP){
        rant[d] <- upS[d,]%*%cbind(upccin[d,])
       }
         marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])-rant
         
      if (length(rows)>1){
         sumt2<-crossprod(marki)-marki[length(marki)]^2
         that<-0
       for (i in 2: length(marki)){
           that<-that+marki[i-1]*marki[i]
       }
         Ymtotal <- Ymtotal + sumt2
         Ymtotal2 <-Ymtotal2+ that
    }
    }
    return (c(Ymtotal, Ymtotal2))
}
PhiTermARp <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W, Ti=id2, CC=ccin, newbeta=beta[g,], lag=p){
      Ymtotal<-matrix(0, nrow=lag, ncol=lag)
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- PX[rows,]
      upW <- PW[rows,]
      TT<-Ti[rows]
      upccin <- CC[TT]
      marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])-upccin
       if (length(rows)>lag){
      for(q in (lag+1):length(marki)){
        trans <- c()
        for (w in 1:lag){
          trans[w] <-marki[q-w]
        }
        Ymtotal<-Ymtotal+cbind(trans)%*%rbind(trans)
       Ymtotal2<-Ymtotal2+cbind(trans)*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}

PhiTermARpT <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W, PS=S, Ti=id2, CC=ccin, newbeta=beta[g,], lag=p){
      Ymtotal<-matrix(0, nrow=lag, ncol=lag)
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- as.matrix(PX[rows,])
      upW <- as.matrix(PW[rows,])
      upS <- as.matrix(PS[rows,])
      TT<-Ti[rows]
      upccin <- as.matrix(CC[TT, ])
      JP <- length(rows)
      rant <- c()
      for (d in 1: JP){
        rant[d] <- upS[d,]%*%cbind(upccin[d,])
      }
       marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])-rant
       if (length(rows)>lag){
      for(q in (lag+1):length(marki)){
        trans <- c()
        for (w in 1:lag){
          trans[w] <-marki[q-w]
        }
        Ymtotal<-Ymtotal+cbind(trans)%*%rbind(trans)
        Ymtotal2<-Ymtotal2+cbind(trans)*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}

PhiTermARp1 <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W,  newbeta=beta[g,], lag=p){
      Ymtotal<-matrix(0, nrow=lag, ncol=lag)
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- as.matrix(PX[rows,])
      upW <- as.matrix(PW[rows,])
      marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])
    #  cat("lag ", lag, "\n", sep=" ")
       if (length(rows)>lag){
      for(q in (lag+1):length(marki)){
        trans <- c()
        for (w in 1:lag){
          trans[w] <-marki[q-w]
        }
        Ymtotal<-Ymtotal+cbind(trans)%*%rbind(trans)
       Ymtotal2<-Ymtotal2+cbind(trans)*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}

PhiTermARp2 <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W, PS=S, Ti=id2, CC=timerandom, newbeta=beta[g,], lag=p){
  
      Ymtotal<-matrix(0, nrow=lag, ncol=lag)
      Ymtotal2<-0
       for (n in 1: NN){
        rows<-Entry[which(unit==n)]
        TL <- length(rows)
        upZ <- yystar[rows]
        upX <- PX[rows,]
        upW <- PW[rows,]
        upS <- PS[rows,]
        TT<-Ti[rows]
        upccin <- timerandom[TT,]
         Rant <- c()
         for (l in 1:TL){
          Rant[l] <- upS[l,]%*%cbind(upccin[l,])
        }
        marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])-Rant
       if (length(rows)>lag){
      for(q in (lag+1):length(marki)){
        trans <- c()
        for (w in 1:lag){
          trans[w] <-marki[q-w]
        }
        Ymtotal<-Ymtotal+cbind(trans)%*%rbind(trans)
       Ymtotal2<-Ymtotal2+cbind(trans)*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}

PhiTermARp1 <- function(Entry=entry, unit=id1, NN=GU, yystar=upystar,newb=upbbn, PX=X, PW=W,  newbeta=beta[g,], lag=p){
      Ymtotal<-matrix(0, nrow=lag, ncol=lag)
      Ymtotal2<-0
       for (n in 1: NN){
      rows<-Entry[which(unit==n)]
      upZ <- yystar[rows]
      upX <- PX[rows,]
      upW <- PW[rows,]
      marki <- upZ-upX%*%cbind(newbeta)-upW%*%cbind(newb[n,])
    #  cat("lag ", lag, "\n", sep=" ")
       if (length(rows)>lag){
      for(q in (lag+1):length(marki)){
        trans <- c()
        for (w in 1:lag){
          trans[w] <-marki[q-w]
        }
        Ymtotal<-Ymtotal+cbind(trans)%*%rbind(trans)
       Ymtotal2<-Ymtotal2+cbind(trans)*marki[q]
       }
    }
    }
    return (cbind(Ymtotal, Ymtotal2))
}


PhiProp <- function(Ymt=PhiSig, Ymt2=PhiMean,scale=1, Oldphi=phig){
      
        HATPhi <-pseudoinverse(Ymt)*scale
        HATphi<-HATPhi %*%Ymt2
         cat("HATphi ", HATphi, "\n", sep=" ")
       phi2par <- BiNorm2(mu=HATphi, y1=Oldphi[1],Sigma=HATPhi)
       pmean2 <- phi2par[1]
       pvar2 <- sqrt(phi2par[2])
       if (pmean2>(1+30*pvar2)|pmean2<(-1-7*pvar2)){
       phi2.draw <-runif(1, max=0.99, min=-0.99)
       }else{
       phi2.draw <- rtnorm(1, pmean2, pvar2, upper=0.99, lower=-0.99)
       }
       phi1par <- BiNorm1(mu=HATphi, y2=phi2.draw,Sigma=HATPhi)
       pmean1 <- phi1par[1]
       pvar1 <- sqrt(phi1par[2])
       upp <- 0.99-phi2.draw
       loww <- phi2.draw-0.99
        if (pmean1>(upp+35*pvar2)|pmean2<(-loww-8*pvar2)){
       phi1.draw <-runif(1, max=upp, min=loww)
       }else{
       phi1.draw <- rtnorm(1, pmean1, pvar1, upper=upp, lower=loww)
       }
      phi.draw<-rbind(phi1.draw, phi2.draw)
      phi.draw
}

PhiProp12 <- function(Ymt=PhiSig, Ymt2=PhiMean,scale=1, Oldphi=phig, lag=2){
      
        HATPhi <-pseudoinverse(Ymt)*scale
        HATphi<-HATPhi %*%Ymt2
       phi2par <- BiNorm2(mu=HATphi, y1=Oldphi[1],Sigma=HATPhi)
       pmean2 <- phi2par[1]
       pvar2 <- sqrt(phi2par[2])
       if (pmean2>(1+30*pvar2)|pmean2<(-1-7*pvar2)){
       phi2.draw <-runif(1, max=0.99, min=-0.99)
       }else{
       phi2.draw <- rtnorm(1, pmean2, pvar2, upper=0.99, lower=-0.99)
       }
         if (lag==1){ phi2.draw <- 0}
       phi1par <- BiNorm1(mu=HATphi, y2=phi2.draw,Sigma=HATPhi)
       pmean1 <- phi1par[1]
       pvar1 <- sqrt(phi1par[2])
       upp <- 0.99-phi2.draw
       loww <- phi2.draw-0.99
        if (pmean1>(upp+35*pvar2)|pmean2<(-loww-8*pvar2)){
       phi1.draw <-runif(1, max=upp, min=loww)
       }else{
       phi1.draw <- rtnorm(1, pmean1, pvar1, upper=upp, lower=loww)
       }
      phi.draw<-rbind(phi1.draw, phi2.draw)
      phi.draw
}

PhiPropARp <- function(Ymt=PhiSig, Ymt2=PhiMean, scale=1, Oldphi=phig){

         HATPhi <- pseudoinverse(Ymt)*scale
         HATphi <- HATPhi %*%Ymt2
         unitcir=2
         while(any(unitcir>1)){
      phi.draw <- mvrnorm(1, mu=HATphi, Sigma=HATPhi)
      #     cat("phi.draw ", phi.draw, "\n", sep=" ")
         ord=length(phi.draw)
         GG<-matrix(NA, ncol=ord, nrow=ord)
         GG[2:ord, 1:(ord-1)] <- diag(ord-1)
         GG[1,] <- phi.draw
         GG[-1, ord] <- 0
         eigGG <- eigen(GG)$values
         unitcir <- abs(eigGG)
}
         phi.draw
}



# Function: Arrange data for random effects
random<-function(Var, Pred, index){
   J<-length(unique(na.omit(index)))
     RanFix <- matrix(0, nrow=1, ncol=ncol(Var)*ncol(Pred))
  for (j in 1:J){
    W <- as.matrix(Var[which(index==j),])
    a <- kronecker(Diagonal(ncol(Var)), t(as.matrix(Pred[j,])))
    RanFix<-rbind(RanFix, as.matrix(W%*%a))
   }
   TTT<-RanFix[-1,]
}

randomt <- function(Var, Pred, index1, index2, entry){
   index <- na.omit(index1)
    index2 <- na.omit(index2)
   Var <- na.omit(Var)
   Pred <- na.omit(Pred)
   J1<-length(unique(na.omit(index)))
    WWW <- matrix(NA, nrow=length(entry), ncol=ncol(Var)*ncol(Pred))
   for (j in 1:J1){
     rows <- entry[which(index==j)]
     TT <- index2[rows]
     W <- as.matrix(Var[rows,])
     PW <- Pred[TT,]
     Ee <- length(rows)
     Tt <- matrix(NA, ncol=ncol(Var)*ncol(Pred), nrow=Ee)
     for (s in 1:Ee){
       Tt[s,] <- as.vector(t(cbind(W[s,])%*%rbind(PW[s,])))
     }
     WWW[rows,] <- Tt
       
   }
   WWW
}
  

Random.Coef<- function(Cova=UnitCova, random=random.pred, coef.pred=beta.pred, RBB=RB ){
    random <- as.matrix(random)
    coef.pred <- as.matrix(coef.pred)
    N.unit <- nrow(Cova)
    N.coef <- RBB
   beta.random <-matrix(NA, ncol=N.coef, nrow=N.unit)
     for (n in 1: N.unit){
       beta1 <- c()
       for (j in 1:N.coef){
         beta1[j] <- Cova[n,]%*%cbind(coef.pred[j,])+random[n,j]
    }
            beta.random[n,] <- beta1
     }
     return(as.vector(t(beta.random)))
}

Random.Coef.gamma<- function(Cova=UnitCova, random=random.pred, coef.pred=beta.pred, RBB=RB ){
    random <- as.matrix(random)
    coef.pred <- as.matrix(coef.pred)
    N.unit <- nrow(Cova)
    N.coef <- RBB
   beta.random <-c()
     for (n in 1: N.unit){
       beta.random[n] <-   Cova[n,]%*%cbind(coef.pred)+random[n]
    }
     return(beta.random)
}
# To center the data
Rescale <- function(XXX){
  k <- ncol(XXX)
  n <- nrow(XXX)
  X.rescaled<-matrix(NA, ncol=k, nrow=n)
 for (i in 1: k){
  X.center<-XXX[,i]-mean(XXX[,i])
  nw<-sum(X.center^2)/n                
  X.rescaled[,i]<-X.center/sqrt(nw)
}
  X.rescaled
}
# SGT=GT; Sid2=id2; Sid1=id1; Sttn=ttn; nystar=ystar.draw; SNXstar=XXstar; SNWstar=WWstar; Sbetadrawn=betadrawn; SNSig=NSig; Sphigg=phigg; Sc.draw=c.draw; Entry=entry

# SGT=GT; Sid2=id2; Entry=entry; Sid1=id1; Sttn=ttn; nystar=ystar.draw; SNXstar=XXstar; SNWstar=WWstar; Sbetadrawn=betadrawn; SNSig=NSig; Sphigg=phigg; Tc.draw=c.draw; Cprior=draw.hc
### Problem here
SimRandTime <- function(SGT=GT, Sid2=id2, Entry=entry, Sid1=id1, Sttn=ttn, nystar=ystar.draw, SNXstar=XXstar, SNWstar=WWstar, Sbetadrawn=betadrawn, SNSig=NSig, Sphigg=phigg, Tc.draw=c.draw, Cprior=draw.hc){
     Sc.draw <- Tc.draw
  for (j in 1: SGT){
             mu<-c()
             sig<-c()

            row1<-Entry[which(Sid2==j)]  # Those entries with t=j
            H <- length(row1)           # The observations at t period
            iden1<-Sid1[row1]            # The units at t period
            tttn <- Sttn[iden1,]         # The random unit-specific effects of those units

            ### Get the data at time period j with the units included
             yt1<-nystar[ row1]
             Xt1<-SNXstar[row1,]
             Wt1<-SNWstar[row1,]

          ### Matrix including data and unit specific parameters and indicators
          #     TData1<- cbind(iden1, row1, yt1, Xt1, Wt1, tttn)

           ### vector of parameters
          #     param <- c(betadrawn, Sphigg, Sc.draw)
#     cat("J: ",j, "\n", sep=" ")      
       if (j< (SGT-1)){

  
        ### Get the data at time period j+1 with the units included
               yt2<-nystar[row1+1 ]
               Xt2<-SNXstar[row1+1,]
               Wt2<-SNWstar[row1+1,]
          #     TData2 <- cbind(yt2, Xt2, Wt2)
         ### Get the data at time period j+1 with the units included
             yt3<-nystar[ row1+2]
             Xt3<-SNXstar[row1+2,]
             Wt3<-SNWstar[row1+2,]
          #   TData3 <- cbind(yt3, Xt3, Wt3)
               
         if (j==1){

   ##### Scheme A ####
        for (h in 1:H){
          
          ddraw <- TRFirst(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg)
          
           sig[h] <- ddraw[1]
           mu[h]<-ddraw[2] # it shoud be sig[h]*Sphigg[2]*mu22+var1^(-1)*mu11 as a mean, but since when add all row1 units together, mu has to be timed with 1/sig[h], the two operations cancelled each other. But this probably is not correct
         }
         svar<-1/(sum(sig)+Cprior)
        mmu <- svar*sum(mu)
        Sc.draw [j]<-rnorm(1, mmu, sqrt(svar))

      #    ytn1 <- yt1; yt1 <- yt2;  yt2 <- yt3
       #   Xt1<-Xt2; Xt2<-Xt3;  Wt1<-Wt2;  Wt2<-Wt3
          # row0<-row1
         iden0<-iden1
          }

         if (j==2){
            inter2 <- setdiff(iden1, iden0)
             JJ <- c()
             for (i in 1: length(inter2)) JJ[i] <-which(iden1==inter2[i] )
            for (h in JJ){
              
           ddraw2 <- TRFirst(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg)
          #   cat("ddraw2 ",ddraw2, "\n", sep=" ")
           sig[h] <- ddraw2[1]
           mu[h]<-ddraw2[2]  # it shoud be sig[h]*Sphigg[2]*mu22+var1^(-1)*mu11 as a mean, but since when add all row1 units together, mu has to be timed with 1/sig[h], the two operations cancelled each other. But this probably is not correct
         }
               
             inter <- intersect(iden1, iden0)
             FF <- c()
             for (i in 1: length(inter)) FF[i] <-which(iden1==inter[i] )

         for (h in FF){
               row1[which(row1<2)] <-2 
               ytn1<-nystar[row1-1 ]
               Xtn1<-SNXstar[row1-1,]
               Wtn1<-SNWstar[row1-1,]

           ddraw3 <-TRSecond(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg, Cynt1=ytn1, CXnt1=Xtn1, CWnt1=Wtn1)
      
           sig[h] <-ddraw3[1]
           mu[h]<-ddraw3[2]
         }
              
         svar<-1/(sum(sig)+Cprior)
        mmu <- svar*sum(mu)
        Sc.draw [j]<-rnorm(1, mmu, sqrt(svar))
         iden01 <- iden0
         iden0<-iden1
             }

               if (j!=1 & j!=2){
           inter <- setdiff(iden1, iden0)
             FF <- c()
             for (i in 1: length(inter)) FF[i] <-which(iden1==inter[i] )
         ##### Scheme A ####
        for (h in FF){
         ddraw4 <- TRFirst(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg)
          
           sig[h] <- ddraw4[1]
           mu[h]<-ddraw4[2] 
         }
            inter1 <- setdiff(intersect(iden1, iden0), iden01)
             FF1 <- c()
             for (i in 1: length(inter1)) FF1[i] <-which(iden1==inter1[i] )
               row1[which(row1<2)]=2
               ytn1<-nystar[row1-1 ]
               Xtn1<-SNXstar[row1-1,]
               Wtn1<-SNWstar[row1-1,]
              for (h in FF1){
                
                ddraw5 <- TRSecond(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg, Cynt1=ytn1, CXnt1=Xtn1, CWnt1=Wtn1)
      
           sig[h] <-ddraw5[1]
           mu[h]<-ddraw5[2]
         }
            inter2 <- intersect(iden1, iden0)
            inter3 <- intersect(inter2, iden01)
             FF2 <- c()
             for (i in 1: length(inter3)) FF2[i] <-which(iden1==inter3[i] )
                 for (h in FF2){
                   row1 [which(row1<=2)]<-3
               ytn2<-nystar[row1-2 ]
               ytn1<-nystar[row1-1 ]
         
           ddraw6 <- TRThird(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg, Cynt1=ytn1, Cynt2=ytn2)

           sig[h] <-ddraw6[1]
           mu[h]<-ddraw6[2]
             }
          svar<-1/(sum(sig)+Cprior)
        mmu <- svar*sum(mu)
        Sc.draw [j]<-rnorm(1, mmu, sqrt(svar))
         iden01 <- iden0
         iden0<-iden1
           
         }
             }
             
           
    if (j==(SGT-1)){
       
               yt2<-nystar[row1+1 ]
               Xt2<-SNXstar[row1+1,]
               Wt2<-SNWstar[row1+1,]
               
           inter <- setdiff(iden1, iden0)
             FF <- c()
             for (i in 1: length(inter)) FF[i] <-which(iden1==inter[i] )
         ##### Scheme A ####
        for (h in FF){
        ddraw7 <- TRSecondLast(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, CWt1=Wt1, CWt2=Wt2, CXt1=Xt1, CXt2=Xt2, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg)
            sig[h] <-ddraw7[1]
            mu[h]<-ddraw7[2]
         }
               
            inter1 <- setdiff(intersect(iden1, iden0), iden01)
             FF1 <- c()
             for (i in 1: length(inter1)) FF1[i] <-which(iden1==inter1[i] )
               row1[which(row1<2)]=2
               ytn1<-nystar[row1-1 ]
               Xtn1<-SNXstar[row1-1,]
               Wtn1<-SNWstar[row1-1,]
              for (h in FF1){
         ddraw8 <- TRSecondL2(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2, CWt1=Wt1, CWt2=Wt2, CXt1=Xt1, CXt2=Xt2, Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg, Cynt1=ytn1, CXnt1=Xtn1, CWnt1=Wtn1)
         
            sig[h] <-ddraw8[1]
            mu[h]<-ddraw8[2]
         }
            inter2 <- intersect(iden1, iden0)
            inter3 <- intersect(inter2, iden01)
             FF2 <- c()
             for (i in 1: length(inter3)) FF2[i] <-which(iden1==inter3[i] )
                 for (h in FF2){
               row1 [which(row1<=2)]<-3
               ytn2<-nystar[row1-2 ]
               Xtn2<-SNXstar[row1-2,]
               Wtn2<-SNWstar[row1-2,]
                 
               ytn1<-nystar[row1-1 ]
               Xtn1<-SNXstar[row1-1,]
               Wtn1<-SNWstar[row1-1,]
         
           ddraw9 <-TRSecondL3(TT=j, c.drawn=Sc.draw, Cyt1=yt1, Cyt2=yt2,  CWt1=Wt1, CWt2=Wt2,  CXt1=Xt1, CXt2=Xt2,  Cbeta=Sbetadrawn, Cb=tttn, CNSig=SNSig, unit=h, Cphi=Sphigg, Cynt1=ytn1, Cynt2=ytn2)

            sig[h] <-ddraw9[1]
            mu[h]<-ddraw9[2]
             }
         svar<-1/(sum(sig)+Cprior)
        mmu <- svar*sum(mu)
        Sc.draw [j]<-rnorm(1, mmu, sqrt(svar))
         iden01 <- iden0
         iden0<-iden1
           
             }

           if (j==SGT){
           inter <- setdiff(iden1, iden0)
             FF <- c()
             for (i in 1: length(inter)) FF[i] <-which(iden1==inter[i] )
         ##### Scheme A ####
        for (h in FF){
           sig[h]<-SNSig[1,1]^(-1)
         mu[h]<-sig[h]*(yt1[h]-Xt1[h,]%*%cbind(Sbetadrawn)-Wt1[h,]%*%cbind(tttn[h,]))
       }
            inter1 <- setdiff(intersect(iden1, iden0), iden01)
             FF1 <- c()
             for (i in 1: length(inter1)) FF1[i] <-which(iden1==inter1[i] )
               row1[which(row1<2)]=2
               ytn1<-nystar[row1-1 ]
               Xtn1<-SNXstar[row1-1,]
               Wtn1<-SNWstar[row1-1,]
              for (h in FF1){
          cmu1<-ytn1[h]-Xtn1[h,]%*%cbind(Sbetadrawn)-Wtn1[h,]%*%cbind(tttn[h,])
         cmu2<-ytn1[h]-Xtn1[h,]%*%cbind(Sbetadrawn)-Wtn1[h,]%*%cbind(tttn[h,])
         mu12<-c(cmu1, cmu2)
          x1norm<-BiNorm2(mu=mu12, Sigma=SNSig, y1=Sc.draw[1])
          mu11<-x1norm[1]
          var1<-x1norm[2]

          sig[h] <- var1^(-1)
           mu[h]<-var1^(-1)*mu11
         }
            inter2 <- intersect(iden1, iden0)
            inter3 <- intersect(inter2, iden01)
             FF2 <- c()
             for (i in 1: length(inter3)) FF2[i] <-which(iden1==inter3[i] )
                 for (h in FF2){
               row1 [which(row1<=2)]<-3
               ytn2<-nystar[row1-2 ]
               ytn1<-nystar[row1-1 ]
         
            mu21<-(yt1[h]-Sphigg[1]*ytn1[h]-Sphigg[2]*ytn2[h])- Xt1[h,] %*% cbind(Sbetadrawn)- Wt1[h,]%*%cbind(tttn[h,])+Sphigg[1]*Sc.draw[j-1]+Sphigg[2]*Sc.draw[j-2]
           sig[h] <-1
           mu[h]<-mu21
             }
       svar<-1/(sum(sig)+Cprior)
        mmu <- svar*sum(mu)
        Sc.draw [j]<-rnorm(1, mmu, sqrt(svar))
         iden01 <- iden0
         iden0<-iden1
         }
 # cat("sig: ",sig, "\n", sep=" ")
 # cat("mu: ", mu, "\n", sep=" ")
}
    Sc.draw
}

# Function: to compute the value of a multivariate normal density kernel 
 ss<-function(ff, Yb, Xb, bb){
     dds<--1/2*(t(Yb-Xb%*%cbind(bb))%*%pseudoinverse(ff)%*%(Yb-Xb%*%cbind(bb))) 
        dds<-as.vector(dds)
     }

# function: sum the term for all the observations --- for use of posterior variance (normal) update
sumtrans<-function(Mat){
         total<-0
          for (r in 1: nrow(Mat)){
               total<-total+outerself(Mat[r,])
          }
         total
       }

# Another transformation for sampling of \rho
trans2<-function(YY, XX, BB){
  YYstar<-c()
      for (i in 1: length(YY)){
        YYstar[i]<-YY[i]-XX[i, ]%*%BB
      }
  YYstar
 }



##### This is the function for the first time period but not the second last of each unit appearing in the dataset. 
# TT=j; c.drawn=c.draw; Cyt1=yt1; Cyt2=yt2; Cyt3=yt3; CWt1=Wt1; CWt2=Wt2; CWt3=Wt3; CXt1=Xt1; CXt2=Xt2;  CXt3=Xt3; Cbeta=betadrawn; Cb=tttn; CNSig=NSig; unit=h; Cphi=phigg

TRFirst <- function(TT=j, c.drawn=c.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=betadrawn, Cb=tttn, CNSig=NSig, unit=h, Cphi=phigg){
         cmu1<-Cyt1[unit]-CXt1[unit,]%*%cbind(Cbeta)-CWt1[unit,]%*%cbind(Cb[unit,])
         cmu2<-Cyt2[unit]-CXt2[unit,]%*%cbind(Cbeta)-CWt2[unit,]%*%cbind(Cb[unit,])
         mu12<-c(cmu1, cmu2)
         
          x1norm<-BiNorm1(mu=mu12, Sigma=CNSig, y2=c.drawn[2])
          mu11<-x1norm[1]
          var1<-x1norm[2]
         
          mu22<-c.drawn[TT+2]-Cphi[1]*c.drawn[TT+1]-(Cyt3[unit]-Cphi[1]*Cyt2[unit]-Cphi[2]*Cyt1[unit])+CXt3[unit,]%*% cbind(Cbeta)+ CWt3[unit,]%*%cbind(Cb[unit,])
          cvar <- var1^(-1)+Cphi[2]^2
          cmean <- var1^(-1)*mu11+Cphi[2]*mu22
          return(c(cvar, cmean))
         }



TRSecond <- function(TT=j, c.drawn=c.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=betadrawn, Cb=tttn, CNSig=NSig, unit=h, Cphi=phigg, Cynt1=ytn1, CXnt1=Xtn1, CWnt1=Wtn1){
  
     cmu1<-Cynt1[unit]-CXnt1[unit,]%*%cbind(Cbeta)-CWnt1[unit,]%*%cbind(Cb[unit,])
         cmu2<-Cyt1[unit]-CXt1[unit,]%*%cbind(Cbeta)-CWt1[unit,]%*%cbind(Cb[unit,])
         mu12<-c(cmu1, cmu2)
          x1norm<-BiNorm2(mu=mu12, Sigma=CNSig, y1=c.drawn[1])
          mu11<-x1norm[1]
          var1<-x1norm[2]
       
          mu22<-c.drawn[TT+1]-Cphi[2]*c.drawn[TT-1]-(Cyt2[unit]-Cphi[1]*Cyt1[unit]-Cphi[2]*Cynt1[unit])+CXt2[unit,]%*% cbind(Cbeta)+ CWt2[unit,]%*%cbind(Cb[unit,])
      
          mu23<-c.drawn[TT+2]-Cphi[1]*c.drawn[TT+1]-(Cyt3[unit]-Cphi[1]*Cyt2[unit]-Cphi[2]*Cyt1[unit])+CXt3[unit,]%*% cbind(Cbeta)+ CWt3[unit,]%*%cbind(Cb[unit,])
        
         cvar <- var1^(-1)+Cphi[2]^2+Cphi[1]^2
          cmean <- var1^(-1)*mu11+Cphi[2]*mu23+Cphi[1]*mu22
          return(c(cvar, cmean))
         }

TRSecondLast <- function(TT=j, c.drawn=c.draw, Cyt1=yt1, Cyt2=yt2, CWt1=Wt1, CWt2=Wt2, CXt1=Xt1, CXt2=Xt2, Cbeta=betadrawn, Cb=tttn, CNSig=NSig, unit=h, Cphi=phigg){
         cmu1<-Cyt1[unit]-CXt1[unit,]%*%cbind(Cbeta)-CWt1[unit,]%*%cbind(Cb[unit,])
         cmu2<-Cyt2[unit]-CXt2[unit,]%*%cbind(Cbeta)-CWt2[unit,]%*%cbind(Cb[unit,])
         mu12<-c(cmu1, cmu2)
         
          x1norm<-BiNorm1(mu=mu12, Sigma=CNSig, y2=c.drawn[2])
          mu11<-x1norm[1]
          var1<-x1norm[2]
          cvar <- var1^(-1)
          cmean <- var1^(-1)*mu11
          return(c(cvar, cmean))
         }

TRSecondL2 <- function(TT=j, c.drawn=c.draw, Cyt1=yt1, Cyt2=yt2, CWt1=Wt1, CWt2=Wt2, CXt1=Xt1, CXt2=Xt2, Cbeta=betadrawn, Cb=tttn, CNSig=NSig, unit=h, Cphi=phigg, Cynt1=ytn1, CXnt1=Xtn1, CWnt1=Wtn1){
  
        cmu1<-Cynt1[unit]-CXnt1[unit,]%*%cbind(Cbeta)-CWnt1[unit,]%*%cbind(Cb[unit,])
         cmu2<-Cyt1[unit]-CXt1[unit,]%*%cbind(Cbeta)-CWt1[unit,]%*%cbind(Cb[unit,])
         mu12<-c(cmu1, cmu2)
          x1norm<-BiNorm2(mu=mu12, Sigma=CNSig, y1=c.drawn[1])
          mu11<-x1norm[1]
          var1<-x1norm[2]
       
          mu22<-c.drawn[TT+1]-Cphi[2]*c.drawn[TT-1]-(Cyt2[unit]-Cphi[1]*Cyt1[unit]-Cphi[2]*Cynt1[unit])+CXt2[unit,]%*% cbind(Cbeta)+ CWt2[unit,]%*%cbind(Cb[unit,])
        
         cvar <- var1^(-1)+Cphi[1]^2
          cmean <- var1^(-1)*mu11+Cphi[1]*mu22
          return(c(cvar, cmean))
         }

TRSecondL3 <- function(TT=j, c.drawn=c.draw, Cyt1=yt1, Cyt2=yt2,  CWt1=Wt1, CWt2=Wt2,  CXt1=Xt1, CXt2=Xt2,  Cbeta=betadrawn, Cb=tttn, CNSig=NSig, unit=h, Cphi=phigg, Cynt1=ytn1, Cynt2=ytn2){
  
      mu21<-(Cyt1[unit]-Cphi[1]*Cynt1[unit]-Cphi[2]*Cynt2[unit])- CXt1[unit,] %*% cbind(Cbeta)- CWt1[unit,]%*%cbind(Cb[unit,])+Cphi[1]*c.drawn[TT-1]+Cphi[2]*c.drawn[TT-2]

          mu22<-c.drawn[TT+1]-Cphi[1]*c.drawn[TT-1]-(Cyt2[unit]-Cphi[2]*Cyt1[unit]-Cphi[2]*Cynt1[unit])+ CXt2[unit,]%*% cbind(Cbeta)+CWt2[unit,] %*%cbind(Cb[unit,])

         cvar <- 1+Cphi[1]^2
          cmean <- mu21+Cphi[1]*mu22
          return(c(cvar, cmean))        
 
         }


# TT=j; c.drawn=c.draw; Cyt1=yt1; Cyt2=yt2; Cyt3=yt3; CWt1=Wt1; CWt2=Wt2; CWt3=Wt3; CXt1=Xt1; CXt2=Xt2;  CXt3=Xt3; Cbeta=betadrawn; Cb=tttn; CNSig=NSig; unit=h; Cphi=phigg; Cynt1=ytn1; Cynt2=ytn2
TRThird <- function(TT=j, c.drawn=c.draw, Cyt1=yt1, Cyt2=yt2, Cyt3=yt3, CWt1=Wt1, CWt2=Wt2, CWt3=Wt3, CXt1=Xt1, CXt2=Xt2,  CXt3=Xt3, Cbeta=betadrawn, Cb=tttn, CNSig=NSig, unit=h, Cphi=phigg, Cynt1=ytn1, Cynt2=ytn2){
  
      mu21<-(Cyt1[unit]-Cphi[1]*Cynt1[unit]-Cphi[2]*Cynt2[unit])- CXt1[unit,] %*% cbind(Cbeta)- CWt1[unit,]%*%cbind(Cb[unit,])+Cphi[1]*c.drawn[TT-1]+Cphi[2]*c.drawn[TT-2]

          mu22<-c.drawn[TT+1]-Cphi[1]*c.drawn[TT-1]-(Cyt2[unit]-Cphi[2]*Cyt1[unit]-Cphi[2]*Cynt1[unit])+ CXt2[unit,]%*% cbind(Cbeta)+CWt2[unit,] %*%cbind(Cb[unit,])

          mu23<-c.drawn[TT+2]-Cphi[1]*c.drawn[TT+1]-(Cyt3[unit]-Cphi[1]*Cyt2[unit]-Cphi[2]*Cyt1[unit])+ CXt3[unit,]%*% cbind(Cbeta)+CWt3[unit,]%*%cbind(Cb[unit,])

         cvar <- 1+Cphi[2]^2+Cphi[1]^2
          cmean <- mu21+Cphi[2]*mu23+Cphi[1]*mu22
          return(c(cvar, cmean))        
 
         }

Wishart <- function ( element, col){
   mm <- array(0, c(col,col))
   mm[lower.tri(mm, diag = TRUE)] <-element
   res <- mm + t(mm)
   diag(res) <- diag(res)/2
   return (res)
 }
#-----------------------------------------------------------------------
V.matrix <- function(Omega=covariance, col=col){
  kap <- min(eddcmp(Omega)$evalues)/2
  J <- Omega-diag(kap, col)
  V <- t(chol(J))
  return (list(kap, V))
}





 VaryingCoef <- function(KK, Betai, Residual, m, RB, GU, unitint.add){
   betaii <- matrix(NA, nrow=1, ncol=RB*GU)
   for (i in 1: m){
     # if there is a random intercept
     if (unitint.add==1){
        Res1 <- matrix(Residual[i,], ncol=(RB+1))
        Res <- as.matrix(Res1[, -1])
     }else{
        Res <- matrix(Residual[i,], ncol=(RB))
     }
     Betaii <- matrix(Betai[i,], nrow=RB)
     beta2 <- Random.Coef(Cova=KK, random=Res, coef.pred=Betaii, RBB=RB)
     betaii <- rbind(betaii, as.vector(beta2))
   }
   betai <- betaii[-1,]
   return(betai)
 }

 MultifactorCoef <- function(Zt, Betai, theta, gamma, m, RB, GU){
    betaii <- matrix(NA, nrow=GT, ncol=GU)
   for (i in 1: m){
     mtheta <- matrix(theta[i,], ncol=RB)
     mgamma <- matrix(gamma[i,], ncol=RB)
     mbetai <- matrix(betai[i,], ncol=RB)
     for (j in 1: GU){
       for (q in 1:GT){
         betaii[q,j] <- Zt[q,]%*%mbetai[j,]+mtheta[q,]%*%mgamma[j,]
       }
     return(betaii)
     }
   }
  }
    
